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PREFACE 


The  research  described  in  this  report  was  funded  by  Grant  No.  MSM 
9224488  from  the  National  Science  Foundation  and  by  Grant  No.  F49620-93-1-0435 
from  the  Air  Force  Office  of  Scientific  Research.  The  report  is  divided  into  two 
parts.  Part  1  contains  the  theoretical  development,  implementation  and  results 
of  a  lattice  type  model.  Part  II  contains  an  experimental  verification  of  the  lattice 
type  model  through  photoelastic  experiments  and  image  processing  of 
photosensitive  disks  tested  in  a  simple  shear  apparatus  developed  as  part  of  this 
research.  Each  Part  is  self  contained. 
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ABSTRACT 


The  objective  of  this  research  is  to  develop  a  lattice  type  model  for 
describing  the  micro  mechanical  behavior  of  particulate  media.  In  this  model, 
each  particle  is  replaced  by  a  lattice  and  the  interactions  between  each  particle 
and  its  neighbors  are  described  through  contact  nodes.  An  assembly  of  particles 
is  then  transformed  into  a  two  dimensional  truss.  Bars  connecting  the  contact 
points  of  the  particles  represent  the  particles  and  constraints  are  prescribed  at  the 
nodes  of  the  truss  to  describe  sliding,  bonding  and  rebonding.  When  a 
deformation  mechanism  develops,  the  truss  becomes  a  pin  jointed  frame  to  derive 
the  kinetics  of  the  assembly. 

The  mechanics,  formulation  and  numerical  implementation  of  the  lattice 
type  model  are  depicted.  Numerical  tests  comprising  two  dimensional  assemblies 
of  disks,  arranged  as  very  loose  and  very  dense  packings,  were  subjected  to 
simple  shear  deformation.  The  results  show  the  distribution  of  die  internal  forces, 
the  load  and  displacement  paths,  and  the  deformation  mechanisms.  Although  the 
packings  were  homogenous  and  isotropic,  the  applied  simple  shear  strains  result 
in  inhomogeneous  deformation.  Deformation  takes  place  through  chains  or 
clusters  of  particles.  Photoelastic  experiments  on  disks,  presented  in  Part  n  of  this 
report,  gave  results  that  are  in  agreement  with  the  predictions  from  the  lattice 
type  model.  The  lattice  type  model  is  conceptually  simple  but  has  some  powerful 
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features  (for  examples,  initial  imperfections,  anisotropy,  and  particle  crushing  can 
be  considered)  that  can  provide  significant  insights  into  the  micro  mechanical 
behavior  of  particulate  media. 
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CHAPTER  1 


INTRODUCTION 


1.1  Particulate  materials 


A  particulate  material  may  be  defined  as  one  composed  of  mutually 
contacting  solid  particles  or  structural  units  within  the  liquid  and/or  the  gaseous 
phase  and  which  exhibit  dilatancy,  contractancy  and  sensitivity  to  applied 
stresses.  The  description  of  a  particle  as  "solid"  or  "grain"  is  relative.  A  'solid 
particle'  is  that  part  of  the  skeleton  which  moves  as  a  single  imit  during  a 


Fig.  1.1  Structure  of  a  Particulate  material 
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mechanical  loading  process.  Thus  a  grain  may  be  defined  as  an  assemblage  of 
elementary  particles.  The  structure  of  a  particulate  material  is  schematically 
shown  in  Fig.  1.1. 

The  grains  of  a  particulate  material  are  typically  in  mutual  contact.  The 
existence  of  the  mutual  contact  restricts  the  freedom  or  autonomy  of  motion  of 
the  individual  grains  and  thus  conditions  the  strength  and  the  rigidity  of  a 
particulate  material.  The  strength  and  rigidity  depend  on  the  number  and 
strength  of  the  contact  bonds  which  are  themselves  a  function  of  the  size,  shape, 
roughness,  relative  position  and  strength  of  the  grains,  and  of  the  nature  of 
interaction  between  the  grains  (Feda,  1982).  The  most  significant  feature  of 
particulate  materials  is  that  their  deformation  is  brought  about  by  mutual  sliding 
of  the  grains  (intergranular  deformation),  in  contrast  to  the  deformation  of  the 
individual  grains  (intragranular  deformation)  typical  of  continuous  media.  Thus 
the  limitations  of  the  continuum  theories  for  the  behavior  of  particulate  materials 
which  do  not  account  for  the  micro-structural  effects  have  led  to  tiie  development 
of  models  dependent  upon  particle  level  interactions. 

1.2  Mechanics  of  Particulate  Materials 

In  a  structural  approach  to  the  mechanics  of  particulate  materials,  the 
grains  are  idealized  to  a  set  of  elastic  spheres  in  mutual  contact.  Most  of  the 


research  on  the  behavior  of  a  set  of  elastic  spheres  is  based  on  Hertz's  solution 


(Johnson,  1985)  for  two  deformable  elastic  bodies  in  contact.  Hertz's  theory  of 
normal  contacts  was  extended  by  Mindlin  (1949)  to  include  the  effect  of  tangential 
forces  at  the  contacts.  The  Hertz-Mindlin  theory  of  contacts  provides  a  foundation 
for  research  into  the  mechanical  behavior  of  particulate  media. 

The  uncertainties  regarding  interior  stresses  in  samples  of  sand  led  to  the 
development  of  models  of  granular  media  that  are  geometrically  simpler  than 
sand.  These  models  consist  of  assemblies  of  disks  (Fig.  1.2)  or  spheres  and  may 
be  analytical,  physical  or  numerical. 


Fig.  1.2  Disk  assembly  model  of  Particulate  medium 


Experimental  difficulties  and  limitations  of  the  analytical  approaches  which 
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will  be  discussed  in  the  next  chapter  have  encouraged  the  development  of 
numerical  techniques.  Cundall  and  Strack  (1979a)  developed  a  powerful 
numerical  technique  called  the  Discrete  Element  Method  (DEM),  wherein  an 
assembly  of  rigid  disks  are  linked  by  normal  and  tangential  springs  at  contacts 
and  the  problem  is  treated  as  a  dynamical  one.  Kishino  (1988)  modified  this 
method  for  the  analysis  of  statical  behavior  of  granular  materials  and  called  it  as 
Disk  Model  Analysis  (DMA).  Expressions  for  elastic  constants  (Bathurst  and 
Rothensburg,  1988),  constitutive  relationships  (Chang  and  Misra,  1989a,b)  and 
structure  of  shear  bands  (Bardet  and  Proubet,  1991b  and  1992)  have  been 
investigated  for  particulate  materials  using  various  concepts /techniques 
incorporating  the  effects  of  micro-structural  properties. 

1.3  Statement  of  the  Problem 

Recently,  there  has  been  a  growing  interest  among  physicists  and  engineers 
(Krajcinovic  and  Basista,  1991;  Adley  and  Sadd,  1992;  Dai  and  Frantziskonis,  1994; 
Frantziskonis  et.al,  1994)  to  study  the  behavior  of  materials  through  a  lattice-t5q?e 
network.  The  mechanical  behavior  of  particulate  materials,  composed  of  granular 
and/or  fibrous  microstructures,  is  inherently  involved  with  the  transmission  of 
loads  along  discrete  paths  within  the  material  (Adley  and  Sadd,  1992).  Based  on 
the  fact  that  granular  materials  transmit  loads  only  through  contacts,  such 
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materials  have  been  modelled  with  network  theories  (Trollope  and  Burman,  1980). 
It  is  this  network  like  pattern  of  load  distribution  in  a  particulate  material  in  the 
form  of  a  lattice  that  has  led  to  the  proposed  development  of  the  lattice  type 
model  for  particulate  media  in  this  research. 

Oda  (1972a,b,c)  studied  the  effects  of  initial  fabrics  and  their  relations  on 
the  mechanical  properties  of  granular  materials.  The  fabric  of  a  granular  material 
is  quantified  by  several  measures  such  as  the  axial  ratios  of  particles,  orientation 
of  contact  normals  and  orientation  of  branch  vectors  (Korushi  and  Naruse,  1988). 
The  lattice  type  model  would  very  well  incorporate  the  effects  of  fabric,  as  the 
lattice  structure  is  evolved  by  linking  the  particle  contacts.  An  attractive  feature 
of  the  lattice  type  model  would  be  the  modelling  of  particles  within  their  micro¬ 
structures  as  different  from  the  modelling  of  contacts  in  the  current  numerical 
approach,  DEM.  Thus  the  lattice  type  model  would  be  a  conceptually  simple 
model  incorporating  several  features  of  particulate  media  such  as  particle  shapes, 
particle  interactions,  particle  kinetics  and  load  path  and  is  expected  to  reveal 
significant  aspects  of  particulate  material  behavior. 

1.4  Objective  of  this  research 

The  objective  of  this  research  is  to  determine  whether  a  lattice  t5q)e 
simulation  can  lead  to  significant  insights  and  imderstanding  of  the  mechanical 
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behavior  of  particulate  media  and  thus  result  in  a  better  understanding  of  it's 

behavior. 

1.5  Scope  of  this  research 

•  Develop  the  mathematical  concepts  of  a  lattice  type  model  for  particulate 
media. 

•  Develop  a  computer  program  for  the  numerical  implementation  of  the 
concepts  formulating  the  lattice  type  model. 

•  Validate  the  potential  of  the  lattice  type  model  by  conducting  numerical 
tests  on  particulate  packings  and  comparing  the  results  with  existing 
laboratory  test  results  on  granular  materials,  DEM  and  experiments  on 
photoelastic  disks. 

•  Explain  the  failure  mechanisms  in  particulate  media  from  numerical  test 
results  and  investigate  the  effects  of  initial  imperfections. 

•  Evaluate  and  interpret  the  physical  quantities  related  to  the  macro¬ 
mechanical  behavior  of  particulate  media  from  the  micro-structural 
properties. 
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1.6  Organization  of  text 

The  introduction  and  background  presented  in  this  chapter  provides  an 
overall  view  of  the  contents  of  this  report.  Chapter.  2  deals  with  the  literature 
review  of  the  research  work  carried  out  to  understand  the  mechanical  behavior 
of  particulate  media.  This  chapter  begins  with  continuum  theories,  a  description 
of  Hertz-Mindlin  contact  theory,  disk  model  analysis  of  particulate  media, 
description  of  discrete  element  method  (DEM),  review  of  the  numerical  tests 
using  DEM,  concept  of  fabric  in  granular  media,  micro-structural  continuum 
theories  and  a  review  of  lattice  modelling  related  to  the  mechanics  of  solids. 

In  Chapter.  3,  the  mathematical  concepts  related  to  the  lattice  type  model 
and  the  numerical  implementation  of  the  concepts  using  a  computer  program  are 
explained.  Chapter.  4  deals  with  numerical  tests  on  particulate  packings  and  the 
discussion  of  the  results.  Chapter.  4  also  includes  a  qualitative  comparison  of  the 
numerical  test  results  with  laboratory  test  results  on  granular  materials, 
evaluation  of  physical  quantities  related  to  macro-mechanical  behavior  of 
particulate  media  and  an  explanation  of  failiire  mechanisms  in  particulate  media. 
The  report  is  concluded  in  Chapter.  5,  listing  out  the  conclusions  of  this  research 
and  recommendations  for  future  research. 
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CHAPTER  2 


LITERATURE  REVIEW 


2.1  Continuum  theories 

Continuum  theories  constitute  the  earliest  attempts  in  the  modeling  of  the 
loading  response  of  particulate /granular  materials.  In  this  macroscopic  theory, 
constitutive  relationships  are  derived  from  the  analysis  of  experimentally 
determined  stress-strain  response  of  granular  materials.  The  prediction  of  loading 
response  requires  (1)  a  theoretical  model  that  exhibits  the  major  features  of  the 
material  response  and  (2)  the  quantitative  determination  of  material  parameters 
which  characterize  the  model. 

The  most  commonly  used  continuum  theories  can  be  broadly  classified 
under  two  categories  as  the  deformation  theory  and  plasticity  theory. 
Deformation  theory  implies  the  existence  of  a  relatively  simple  relationship 
between  the  stress  and  strain  tensors.  The  hyperbolic  soil  model  of  Duncan  and 
Chang  (1970)  is  an  example  of  a  model  which  uses  deformation  theory.  The 
model  is  based  on  a  generalization  of  Hooke’s  law  for  isotropic  behavior  and  it 
cannot  model  anisotropic  responses,  shear  dilatancy  and  postpeak  strain 
softening.  Limitations  of  the  deformation  theory  led  to  the  development  of  a 
number  of  plasticity-based  models.  The  most  widely  used  of  these  models  is  the 
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incremental  theory  of  plasticity  in  which  the  total  strain  increment  is  expressed 
as  the  sum  of  the  elastic  and  plastic  strain  increments. 

A  major  limitation  of  the  continuum  theories  when  applied  to  granular 
materials  is  that  no  accoxmt  is  taken  of  the  effects  of  micro-structural  properties 
and  particle  level  interactions.  This  led  to  the  development  of  models  based  on 
micro-mechanical  behavior  of  particles  for  which  the  basis  is  the  contact  theories 
of  Hertz  (1882)  and  Mindlin  (1949). 


2.2  Hertz-Mindlin  theory  of  contacts 


The  relative  movement  of  two  deformable  elastic  bodies  in  contact  resulting 
from  normal  force  was  studied  by  Hertz  (1882).  When  two  elastic  bodies  are 
pressed  together  with  a  normal  force,  'N'  (Fig.  2.1),  a  small  circular  surface  is 
developed  with  a  radius  called  contact  radius.  Hertz  derived  expressions  for  the 
contact  radius  (a),  the  normal  displacement  (6„)  and  the  normal  stiffriess  (K^)  as 


N 


N 


Fig.  2.1  Hertzian  contact  of  two  elastic  bodies  under  a  normal  load  'N' 
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where  R  is  the  radius  of  the  sphere,  v  is  the  Poisson's  ratio  and  G  is  the  shear 
modulus. 

Mindlin  (1949)  included  the  effect  of  tangential  forces  at  the  contact  and 
found  that  the  tangential  force  at  the  contact  induces  a  deformation  in  the  form 
of  a  slip  over  a  part  of  the  contact  surface.  When  the  tangential  force  exceeds  the 
frictional  strength  at  the  contact,  sliding  takes  place.  Based  on  Mindlin's  (1949) 
results,  the  relative  tangential  displacement  and  the  tangential  contact  stiffness  are 
given  by 
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where  |i  is  the  coefficient  of  friction  at  the  contact,  T  is  the  tangential  force  at  the 
contact,  (6g)  is  the  tangential  displacement  at  the  contact  and  {K^  is  the 
tangential  stiffness  at  the  contact. 

The  expressions  for  normal  stiffness  and  the  tangential  stiffness  derived 
by  Hertz  and  Mindlin  provided  the  basis  for  research  into  the  micro-mechanical 
behavior  of  particulate  media. 

2.3  Disk  model  analysis  of  particulate  media 

The  mechanical  behavior  of  particulate  media,  composed  of  a  random 
packing  of  grains  of  various  sizes,  is  studied  by  geometrically  simpler  models 
consisting  of  assemblies  of  disks  or  spheres.  The  outcome  of  such  research  has 
led  to  a  better  xmderstanding  of  the  behavior  of  particulate  media  and  to  the 
establishment  of  new  concepts  or  models  for  direct  application  in  practice.  The 
solution  technique  for  these  models  could  be  analytical,  physical  or  numerical. 
The  analytical  approach  is  to  evolve  closed  form  solutions.  The  physical  approach 
is  to  conduct  laboratory  experiments.  The  numerical  approach  is  to  conduct 
numerical  tests  through  numerical  computations  based  on  certain  mathematical 
concepts  and  implementing  the  concepts  using  computer  programs. 

An  analytical  model  for  cubic  arrays  of  spheres  of  uniform  size  was 
proposed  by  Deresiewicz  (1958)  which  predicted  a  non-linear  and  hysteretic 
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stress-strain  behavior.  Ultimate  failure  was  also  accommodated  in  the 
formulation.  This  analytic  approach,  however,  was  restricted  to  cubic/hexagonal 
arrays  of  spheres  and  the  spheres  were  uniform  in  size. 

Experiments  on  assemblies  of  optically  sensitive  disks  has  led  to  the  direct 
measurement  of  contact  forces  and  displacements  of  individual  particles.  Forces 
at  contacts  have  been  obtained  from  the  birefringent-effect  in  polarized  light  and 
displacements  were  foimd  by  comparison  of  photographs  made  at  various  stages 
of  the  experiment.  This  method  originates  from  the  work  of  Dantu  (1957)  and 
Wakabayashi  (1957)  and  was  applied  by  De  Josselin  De  Jong  and  Verruijt  (1969), 
Drescher  and  De  Josselin  De  Jong  (1972)  and  Oda  and  Konishi  (1974).  Some 
limitations  of  the  experimental  techniques  are:  (1)  the  material  of  the  disks  should 
be  optically  sensitive  (2)  the  disks  do  not  adequately  represent  real  particulate 
media. 

With  the  advent  of  powerful  computers,  numerical  techniques  have  become 
the  most  widely  used  approach  for  modeling  assemblies  of  disks  and  spheres. 
In  a  numerical  test,  any  data  such  as  stress,  strain,  particle  displacements  and 
packing  structure  at  various  stages  of  the  test  can  be  stored  numerically  and  can 
be  visualized  using  graphics  programs.  The  flexibility  of  numerical  modeling 
extends  to  loading  configurations,  particle  sizes,  size  distributions  and  physical 
properties  of  the  particles. 

Serrano  and  Rodriguez-Ortiz  (1973)  developed  a  numerical  model  for 
assemblies  of  disks  and  spheres.  Hertzian-type  contact  compliances  were  used 
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for  normal  forces,  the  effects  of  tangential  forces  were  considered  according  to  the 
theories  of  Mindlin  and  Deresiewicz  (1953),  and  shape  changes  were  negligible. 
The  numerical  model  was  based  on  matricial  relationships  between  contact  forces, 
particle  displacements  and  the  boimdary  forces  or  displacements.  A  stiffness 
matrix  was  constructed  that  takes  into  account  the  geometrical  arrangement  of 
particles  and  the  current  stiffness  at  each  contact.  Inverting  the  matrix  resulted 
in  the  computation  of  incremental  displacements  from  the  last  known  forces  with 
an  iteration  procedure  to  deal  with  slip  at  contacts.  Only  one  contact  was  allowed 
to  slip  at  a  time  and  it  was  necessary  to  reform  the  stiffness  matrix  when  new 
contacts  were  established  or  contacts  separated.  A  major  limitation  of  the 
approach  was  the  large  computational  time  required  in  solving  the  equations  even 
for  a  very  small  number  (100)  of  particles.  This  led  to  the  development  of 
discrete  element  method  (DEM),  in  which  the  particles  were  treated  as  discrete 
elements  and  the  solution  technique  was  at  particle  level  through  an  iterative 
process. 

2.4  Discrete  Element  Method 

Cundall  and  Strack  (1979a)  developed  a  numerical  technique  called  the 
distinct  element  method  using  a  time  marching  scheme,  based  on  the  idea  that 
time  step  can  be  chosen  small  enough  that  during  a  single  time  step,  disturbances 
cannot  propagate  from  any  disk  further  than  it's  immediate  neighbors.  This 
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method  has  been  successfully  applied  to  the  study  of  granular  media  and  later 
researchers  refer  to  this  method  as  discrete  element  method  (DEM).  The  method 
is  based  on  the  Hertz-Mindlin  contact  theories  and  Newton’s  second  law  of 
motion.  The  elastic  behavior  of  a  disk  is  replaced  by  tire  action  of  springs 
attached  virtually  at  the  contact  points  in  the  normal  and  tangential  directions, 
and  the  disk  itself  is  displaced  and  rotated  rigidly  (Fig.  2.2).  Disks  can  overlap 
each  other  and  only  at  that  time  a  contact  force  is  produced  between  them.  The 
boundary  can  be  either  strain-controlled  or  stress-controlled,  and  the  interaction 
between  a  boundary  element  and  a  disk  is  specified  in  a  similar  way  as  the 
interaction  between  disks.  The  DEM  was  implemented  using  a  computer 
program  called  BALL  and  was  used  to  perform  numerical  tests  on  two- 
dimensional  assemblies  of  disks  under  a  variety  of  boundary  conditions. 


Normal  Tangential 


Fig.  2.2  Modeling  of  disks  in  the  Discrete  Element  Method 
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The  BALL  program  uses  an  explicit  time-marching  scheme  to  track  the 
motion  of  individual  particles  with  the  motion  of  each  being  governed  by 
conservation  of  momentum.  The  sequence  of  calculations  in  the  BALL  program 
is  shown  in  Fig.  2.3.  The  principles  involved  in  the  calculations  can  be  explained 
with  respect  to  two  particles  x  and  y  in  contact  at  point  C  (Fig.  2.3).  Two  unit 
vectors  e,-  and  tj  are  defined  as  pointing  from  the  center  of  disk  x  to  the  center 
of  disk  y  and  rotated  90°  clockwise  respectively.  The  relative  velocities  X),  of 
the  particle  contacts  are  expressed  as  the  components  of  the  linear  and  angular 
velocities  of  the  contacting  particles.  The  normal  and  shear  components  of  these 
velocities  may  then  be  calculated  as  the  projections  of  Xy  in  the  appropriate 
directions  of  6/  and  tj  respectively.  The  velocities  are  assumed  constant  over  a 
time  step /increment  At.  Thus,  the  incremental  displacements  are  obtained  by 
time  integration.  The  normal  and  shear  stiffnesses  and  are  then  used  to 
obtain  the  incremental  normal  and  shear  components  of  the  contact  forces.  These 
are  then  added  to  the  previous  values  to  obtain  total  normal  and  shear  forces 
with  the  constraint  that  the  shear  force  is  limited  to  a  fraction  of  the  normal  force 
(tanOy)  pl^g  a  constant  cohesion.  This  sequence  of  calculations  is  performed  for 
each  contact  prior  to  the  application  of  the  motion  equations  for  any  particle. 
Since  all  forces  are  known  at  every  contact  point  translational  and  rotational 
accelerations  are  calculated  through  a  difference  approximation  to  Newton's 


second  law. 
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n  =  s  =  Xftj 

Relative Displalcements:  An  =  nAt,  As  =  sAt 
Coniaa  Force  Increments:  AF„  =  AF5  =  kjAs 

Total  Forces:  F„  =  F„"'‘ t  AF,,  F,  =  F,""' -r  AF, 
Check  for  Slip:  F,  =  min(Fs,  F„  tancj),,  -f  c) 

Calculate  Moment:  M^=2FjRx,  My=SFyRy 
EQUATIONS  OF  MOTION 


Assume  Force  &  Moment  Constant  over  At:  (f'^  to  t"*''^) 
Acceleration:  Xi“  =  2  Fj/  m  ,  0"  =  2  M  / 1 


3.  Velocity:  At,  (e")At 

4.  Assume  Velocities  constant  over  At:  (t"  to  t""^*) 

5.  Displacements:  X;"*‘  =  X;”  +  (xf'^'^At 

■6.  Rotation:  0""'  =  0"  +  (^’''")At 


y 

rNCRKMENT  TIME' 


Fig.  2.3  Calculation  sequence  in  DEM  (Strack  and  Cundall,  1978) 
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Time  f  ®  refers  to  time  step  t  =  0,  time  t "  refers  to  time  t  =  n  At  and  time 
refers  to  time  t  =  n  At  +  At/2.  The  moments  and  forces  are  assumed  to  be 
constant  from  time  to  so  that  the  velocities  centered  at  time  are 
obtained  by  integrating  the  accelerations  over  At  and  incrementing  the  previous 
values.  A  second  integration  yields  the  new  position  and  rotation  centered  at 
The  time  step  is  taken  to  be  small  enough  so  that  the  particles  are 
unaffected  by  motions  more  than  one  diameter  away.  The  new  total  velocities 
will  then  be  used  to  calculate  the  accelerations  at  time  after  the  force- 
displacement  relations  have  again  been  applied  to  the  contacts.  As  interparticle 
overlaps  change,  contacts  are  deleted  and  new  ones  are  added  and  consequently 
the  changes  in  fabric  occurring  in  granular  media  are  captured.  Cimdall  and 
Strack  (1979b)  modified  the  BALL  program  for  the  modeling  of  particles  in  three 
dimensions  giving  it  a  new  name  called  TRUBAL.  The  computer  program  BALL 
and  the  numerical  tests  conducted  using  it  are  elaborately  described  in  Strack  and 
Cimdall  (1978)  and  Cundall  and  Strack  (1979b).  Several  aspects  of  particulate 
material  behavior  have  been  investigated  using  BALL,  TRUBAL  and  later 
modifications  of  the  original  program  and  concepts  by  many  researchers,  which 
are  briefly  reviewed  in  the  next  section. 

2.5  Review  of  Numerical  Experiments 


Cundall  and  Strack  (1979a)  performed  numerical  biaxial  tests  on  two- 
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dimensional  assemblies  of  disks  using  the  DEM  and  validated  the  potential  of  the 
technique  by  qualitatively  comparing  the  contact  force  vector  plots  (Fig.  2.4)  with 
the  results  of  experiments  on  photoelastic  disks  by  De  Josselin  De  Jong  and 
Verruijt  (1969).  Cundall  and  Strack  (1979c)  performed  numerical  simple  shear  test 
on  a  two-dimensional  assembly  of  cylinders  similar  to  that  performed 
experimentally  by  Oda  and  Konishi  (1974).  Even  though  the  packing 
arrangement  in  the  numerical  test  was  different  from  the  experiment,  the  total 
number  of  particles  and  the  size  distribution  of  the  particles  (ratio  of  number  of 
particles  of  different  sizes)  used  in  the  numerical  test  were  same  as  in  the 
experiment.  The  good  agreement  of  the  plots 


(a)  Ball  program  (b)  Photoelastic  experiment 

Fig.  2.4  Vector  plot  of  contact  forces  from  Ball  program  and 
De  Josselin  De  Jong  (1969) 


of  the  total  shear  force  and  volumetric  strain  versus  shear  strain  between  the 
numerical  test  and  the  experiment  led  to  the  validation  of  DEM  quantitatively. 
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Thus  the  method  was  established  as  a  powerful  tool  for  developing  constitutive 
laws  for  idealized  granular  materials. 

The  mechanisms  observed  in  the  deformation  and  failure  of  granular  media 
was  explained  by  Cimdall  et.  al  (1982).  The  observed  mechanisms  showed  that 
in  granular  media,  forces  are  never  transmitted  across  a  sample  in  a  uniform  way, 
but  are  concentrated  in  'chains'  of  particles.  Between  the  chains,  particles  may 
carry  little  or  no  load.  Sliding  contacts  are  almost  never  observed  between 
particles  comprising  the  major  force  chains.  The  total  number  of  contacts  in  a 
sample  is  foimd  to  decrease  when  the  sample  is  loaded  deviatorically.  The 
contact  deletions  occur  mainly  for  contacts  with  contact  normals  in  the  direction 
of  minor  principal  strain.  The  microscopic  observations  in  the  numerical  tests 
showed  that  a  granular  assembly  consisted  of  two  phases;  Phase  A,  made  up  of 
stiff  columns  which  behaves  like  a  pin-jointed  structure  that  dissipates  no  energy 
and  imposes  local  deformation  patterns  on  Phase  B  particles.  Phase  B  acts  like 
a  truly  plastic  material  that  dissipates  energy  at  many  contacts  and  serves  as  a 
restraint  on  Phase  A. 

Based  on  the  observations  from  the  numerical  tests,  Cxmdall  and  Strack 
(1983)  introduced  a  measure  called  constraint  ratio  (CR)  expressed  in  terms  of  the 
number  of  particles,  number  of  contacts  and  number  of  sliding  contacts.  The 
constraint  ratio  was  used  to  interpret  the  structural  stability  of  granular  materials. 
If  CR  >  1,  the  system  was  formd  to  be  over  constrained  or  statically  indeterminate 
and  the  assembly  was  foimd  to  have  a  positive  resistance  to  deformation.  If  CR 
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<  1,  there  were  more  degrees  of  freedom  than  constraints  and  deformation  was 
found  to  take  place  for  zero  or  negative  increment  in  loading.  If  CR  =  1,  the 
granular  assembly  was  found  to  have  reached  a  condition  of  limiting  stress  and 
at  this  stage  the  assembly  was  foimd  to  develop  a  "mechanism"  and  deformation 
proceeds  in  a  particular  mode  with  no  increase  in  stress. 

Cundall  (1988)  presented  results  from  several  numerical  experiments  on 
dense  sphere  assemblies  in  three  dimensions  using  TRUBAL  and  established 
good  comparison  with  the  results  of  physical  experiments  described  by  Ishibashi 
and  Chen  (1988).  The  nature  of  localization  in  frictional  materials  was 
investigated  by  Cimdall  (1989)  using  two  methods;  (1)  a  continuum  calculation 
with  a  strain-hardening  constitutive  model  and  (2)  using  the  discrete  element 
method.  It  was  shown  that  continuum  calculations  done  with  random  spatial 
distributions  of  physical  properties  can  give  rise  to  localization  in  the  strain¬ 
hardening  regime.  It  was  foimd  that  localized  shear  band  decreases  in  thickness 
as  shearing  proceeds.  The  shear  band  was  characterized  by  dilatation,  reduction 
in  the  parallel  stress  component  and  average  particle  spins  that  deviated 
significantly  from  the  rotation  of  the  bulk  material.  Cundall  et.al  (1989) 
performed  compression/extension  shear  tests  by  applying  strain  rates  to  a  three 
dimensional  granular  assembly  and  evaluated  incremental  elastic  moduli  during 
loading  and  unloading  using  the  discrete  element  method.  This  was  done  with 
an  emphasis  on  the  variation  of  elastic  wave  speeds  as  the  sample  is  loaded  and 
unloaded.  The  elastic  moduli  were  evaluated  in  three  orthogonal  directions  at 
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various  stages  of  loading  and  unloading.  It  was  suggested  that  the  elastic  moduli 
which  is  a  P-wave  measurement  could  be  used  as  a  sensitive  measure  of  strain- 
induced  anisotropy.  The  histograms  of  normal  component  of  contact  forces 
showed  that  force  distributions  were  heavily  biased  towards  the  low  values  and 
it  was  established  that  theories  that  assume  symmetrical  distributions  of  contact 
forces  about  the  mean  were  erroneous. 

Kishino  (1987)  developed  a  numerical  method  called  Disk  Model  Analysis 
(DMA)  for  the  statical  behavior  of  granular  materials.  Unlike  the  DEM  which 
treats  the  problem  as  a  dynamical  one,  the  DMA  uses  a  contact  stiffness  matrix 
for  the  equilibrium  of  the  disks  and  hence  is  more  suitable  for  the  static 
deformation  of  granular  materials.  The  DMA  considers  a  fixing  vector  defined 
as  a  vector  which  consists  of  the  force  and  moment  required  for  fixing  a  disk  at 
it's  center  and  is  calculated  in  terms  of  the  contact  forces  applied  at  it's  contact 
points  and  the  body  force  at  the  center  of  gravity.  The  contact  stiffness  matrix  is 
defined  as  a  set  of  additional  fixing  vectors  required  to  give  unit  movements  of 
a  disk  while  other  disks  were  being  fixed.  To  attain  the  equilibrium  state  of  the 
disk  assembly,  each  disk  is  displaced  and  rotated  iteratively  according  to  the 
contact  stiffness  matrix  defined  from  the  locations  of  neighboring  disks.  The 
DMA  was  used  to  perform  biaxial  tests  on  disk  assemblies.  The  constitutive 
properties  were  evaluated  and  changes  in  fabric  was  evolved.  However,  a 
comparison  of  the  results  between  the  DEM  and  DMA  was  not  presented. 

Thornton  and  Sun  (1993)  performed  numerical  axisymmetric  compression 
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tests  on  dense  and  loose  systems  of  three  dimensional  polydisperse  systems  of 
spheres  using  the  DEM.  The  objective  of  the  study  was  to  evaluate  the  effects  of 
initial  density  and  coefficient  of  interparticle  friction  on  the  evolution  of  various 
micromechanical  parameters.  It  was  observed  that  both  the  maximum  and 
ultimate  angles  of  internal  shearing  resistance  increased  when  the  interparticle 
friction  is  increased.  It  was  noted  that  the  rate  of  dilation  is  significantly  affected 
by  the  magnitude  of  the  interparticle  friction  (p).  The  larger  the  value  of  \i,  the 
higher  is  the  rate  of  dilation,  particularly  for  the  loose  granular  assemblies.  The 
results  of  numerical  simulations  were  compared  with  the  laboratory  axisymmetric 
compression  tests  on  well-rounded  glass  ballotini,  sub-rounded  Welland  river 
sand  and  angular  crushed  glass  performed  by  Parikh  (1967).  The  mobilized  angle 
of  shearing  resistance  (sin  ^  dense  assembly  of  glass  ballotini  with  a  p= 

0.3  was  obtained  as  0.48  and  was  comparable  to  the  sin  obtained  from 
numerical  simulation  as  0.45.  Numerical  test  results  for  dense  assembly  with  a 
p=  0.6  gave  a  sin  =  0.51 .  Results  reported  by  Parikh  (1967)  on  Welland 
river  sand  was  in  the  range  of  0.56  <  sin  reported  by 

Bishop  (1973)  on  Ham  river  sand  was  in  the  range  of  0.54  <  sinO^g^  <  0.68.  The 
numerical  simulation  predicts  a  much  lower  value  for  sin  O^g^^  as  compared  with 
laboratory  test  results  and  it  was  established  that  such  a  high  difference  was  due 
to  the  angularity  of  the  particles  in  the  laboratory  test.  Thus  it  was  concluded  that 
if  the  numerical  simulations  were  to  represent  sand,  it  is  necessary  to  model  non- 
spherical  particles  and  that  merely  increasing  the  coefficient  of  interparticle 
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friction  was  not  sufficient. 

Rothensburg  and  Bathurst  (1992)  presented  results  of  numerical  biaxial 
tests  of  dense  2D  assemblies  of  elliptical  particles.  The  objective  of  the  tests  was 
to  study  the  effects  of  varying  the  aspect  ratio  (A^.  =  larger  dimension/shorter 
dimension  of  ellipse)  of  the  elliptical  particles  and  it  was  shown  that  the  shear 
strength  was  a  function  of  particle  shape.  The  interparticle  friction  was 
maintained  constant  at  a  value  of  p  =  0.5  for  all  the  test  simulations.  For 
1.0  <  1.353,  the  mobilized  angle  of  internal  shear  resistance  for  the 

assembly  increased  with  aspect  ratio  in  the  range  of  26.5°  to  42°.  This  work  also 
demonstrated  the  importance  of  modeling  non-spherical  particles  in  the  computer 
simulations  of  granular  media. 

Bardet  and  Proubet  (1991a)  developed  a  numerical  technique  based  on 
adaptative  dynamic  relaxation  (ADR)  for  the  mechanics  of  granular  materials. 
This  method  was  pursued  as  it  was  more  computationally  efficient  than  the  DEM 
and  it  enhanced  the  rate  of  convergence  and  accuracy.  A  granular  assembly  was 
idealized  as  disks  linked  by  normal  and  tangential  springs  at  contacts  (Fig.  2.2) 
and  the  Hertz-Mindlin  contact  theory  was  assumed  for  the  force-displacement 
law.  The  equilibrium  of  the  disks  was  accomplished  through  the  ADR  technique 
as  different  from  the  Newton's  second  law  in  DEM.  Only  the  numerical  tests 
carried  out  with  the  ADR  technique  is  reviewed  here  and  the  reader  may  refer 
to  Bardet  and  Proubet  (1991a)  for  a  detailed  review  of  the  principles  and 
procedures  involved  in  the  ADR  technique. 
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Bardet  and  Proubet  (1991b)  investigated  the  structure  of  persistent  shear 

bands  in  granular  materials  by  numerically  simulating  an  idealized  two- 

dimensional  assembly  of  particles.  Numerical  biaxial  test  was  performed  on  an 

assembly  of  2000  disks  with  two  rigid  boundaries  in  the  major  principal  stress 

direction  and  two  flexible  boimdaries  in  the  minor  principal  stress  direction.  The 

numerical  test  was  intended  to  simulate  the  rubber  membranes  of  the  triaxial  tests 

in  contrast  to  periodic  boundaries  used  in  the  numerical  test  of  Cundall  (1989). 

The  residual  friction  angle  (0^)  for  the  granular  assembly  was  obtained  as  14° 
with  particles  of  interparticle  friction  angle  of  26.5° .  The  displacement  fields  at 

various  stages  of  axial  strains  indicated  the  formation  of  two  shear  bands  with  an 

inclination  of  52°  and  38°  with  horizontal.  With  the  progress  of  deformation, 

only  the  shear  band  inclined  at  52°  persisted.  This  was  consistent  with  the  Mohr- 

Coulomb  theory  which  predicts  the  failure  plane  at  an  angle  of 

45+0/2  =  45+14/2  =  52°  displacement,  volumetric  strain,  void  ratio, 

rotations  of  particles,  rotations  of  their  neighborhoods  and  contact  orientations 

were  also  examined  inside  the  shear  band.  The  conclusions  of  the  shear  band 

analysis  were  that  shear  band  width  reduces  with  axial  strain  and  the  contacts  are 

oriented  in  the  direction  of  the  shear  band.  The  number  of  these  contacts  per 

particle  is  minimal  and  that  the  rotation  of  particles,  the  gradient  of  their  rotation 

and  the  rotations  of  their  neighborhoods  are  concentrated  inside  the  shear  bands. 

Computer  simulations  using  DEM  and  other  methods  have  been  very 
powerful  techniques  to  understand  the  mechanics  of  granular  materials.  The 
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nvimerical  simulations  which  incorporate  the  microstructural  effects  have  helped 
in  understanding  the  stress-strain  response,  volume  change  behavior, 
development  of  failure  mechanisms,  evaluation  of  physical  quantities  related  to 
the  macro-mechanical  behavior,  stress  and  strain  induced  anisotropy,  formation 
and  persistence  of  shear  bands  and  the  behavior  of  granular  assemblies  in  various 
stress  paths.  A  limitation  of  the  computer  simulations  on  disk/sphere  assemblies 
is  that  they  do  not  incorporate  irregularly  shaped  particles  and  particle  fracture 
and  thus  a  only  limited  capability  exists  to  highlight  and  visualize  specific 
microstructural  events  which  may  control  the  macroscopic  response. 

2.6  Concept  of  Fabric  in  Granular  Media 

Fabric  of  granular  material  means  spatial  arrangement  of  solid  particles 
and  associated  voids  (Oda,  1972a)  or,  in  short,  the  packing  structure.  Oda  (1972a, 
1972b,  1972c)  performed  a  series  of  triaxial  tests  on  granular  materials  of  various 
packing  structures  and  presented  the  concept  of  fabric  in  his  paper  "Significance 
of  fabric  in  granular  mechanics"  (Oda,  1978).  Oda  explained  the  fabric  by  two 
concepts;  (a)  orientation  (b)  packing.  The  orientation  of  the  particles  was  defined 
by  a  term  called  axial  ratio  (Fig.  2.5a)  while  the  packing  was  based  on  the  branch 
vector  and  contact  normal  (Fig.  2.5b). 
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Fig.  2.5  Fabric  measures  of  particles 

Branch  vector  is  the  vector  from  the  centroid  of  particles  to  the  contact 
point.  Contact  normal  is  the  direction  of  the  normal  to  the  contact  surface  when 
two  particles  are  pressed  by  a  contact  force  and  the  direction  of  this  is  obtained 
by  drawing  a  normal  to  the  angle  bisector  of  the  branch  vectors  AC  and  BC  (Fig. 
2.5b).  Quantifying  orientation,  branch  vectors  and  contact  normals  of  the  particles 
in  a  granular  material  in  terms  of  averages,  standard  deviation  and  probability 
density  fimctions  are  elaborately  dealt  in  Oda  (1978).  One  of  the  methods  to 
analyze  the  changes  in  fabric  of  a  granular  material  is  to  construct  a  plot  of  the 
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distribution  of  contact  normals  in  various  directions  in  the  four  quadrants  of  an 
XY  coordinate  system  (Fig.  2.6).  The  shift  in  the  direction  of  the  distribution  of 
contact  normals  indicates  the  shift  in  the  direction  of  the  principal  stress  in  a 
granular  material  as  shown  in  Fig.  2.6. 


(a)  Vertical  load  (b)  Simple  shear 

Fig.  2.6  Frequency  distribution  of  contact  normals  in 
simple  shear  test  (Oda  and  Konishi,  1974) 


The  effect  of  initial  fabric  on  mobilized  strength,  dilatancy  rate  and  secant 
deformation  modulus  were  presented  by  Oda  (1972a)  on  the  results  of  triaxial 
tests.  Based  on  the  experimental  results,  Oda  (1972b  and  1972c)  derived 
expressions  for  the  mobilized  stress  ratio  and  dilatancy  rate  in  terms  of  the  fabric 
measures.  In  order  to  represent  the  distribution  of  contact  normals  as  a 
continuous  function,  Oda  et.al.  (1982)  introduced  a  fabric  tensor  which  is  a 
function  of  the  number  of  contacts  and  the  XY  components  of  the  branch  vectors 


39 


and  the  contact  normals. 

Thornton  and  Barnes  (1986)  performed  biaxial  constant  mean  stress  test 
and  a  constant  volume  test  on  a  dense  assembly  of  1000  disks  using  the  DEM. 
The  objective  of  these  tests  was  to  study  the  induced  structural  anisotropy  during 
the  deformation  of  granular  assemblies.  The  anisotropy  of  a  granular  assembly 
is  defined  by  the  distribution  of  contact  normals  in  the  four  quadrants  of  an  XY 
coordinate  system  (eg.  Fig.  2.6).  If  the  distribution  was  random  and  could  be 
approximated  by  a  imiform  distribution,  then  the  structure  was  in  effect  isotropic. 
Anisotropy  was  indicated  by  a  non-uniform  distribution.  It  was  shown  that  even 
though  the  two  tests  were  subjected  to  different  strain  histories,  the  degree  of 
induced  structural  anisotropy  evolved  in  an  identical  manner.  Also  a  Fourier 
series  expression  was  derived  for  a  probability  density  function  to  express  the 
observed  discrete  distributions  of  contact  normals  (eg.  Fig.  2.6)  as  a  continuous 
distribution.  Nemat-Nasser  and  Mehrabadi  (1983  and  1984)  derived  constitutive 
relationships  based  on  fabric  measures.  As  an  alternative  to  the  fabric  measures 
in  terms  of  the  contact  normals,  Konishi  and  Naruse  (1988)  derived  an  expression 
for  a  term  called  void  tensor  in  terms  of  the  size,  shape  and  orientation  of  local 
voids. 

An  effort  to  verify  the  fabric-stress  relations  through  experiments  was 
carried  out  by  Subhash  et.al  (1991).  Simple  shear  tests  were  conducted  on  a  two- 
dimensional  assembly  of  oval  cross-sectional  rods.  The  components  of  the  fabric 
tensor  were  measured  over  one  cycle  of  shearing.  The  orientations  of  the  principal 
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axes  of  all  commonly  used  fabric  tensors  were  observed  to  change  sharply  with 
the  reversal  of  the  shearing  direction. 

Unlike  the  continuum  models,  the  fabric  based  relationships  incorporated 
the  microstructural  quantities.  Thus  it  is  an  analytical  approach  accounting  for  the 
microstructural  effects  through  statistical  parameters.  Relationships  have  been 
evolved  for  the  various  physical  quantities  related  to  the  macro-mechanical 
behavior  of  particulate  media  in  terms  of  the  initial  fabric  measures.  These 
relationships  require  the  data  consisting  of  particle  geometry,  number  of  contacts 
and  particles  and  contact  orientations  during  various  stages  of  deformation  of  a 
particulate  mediiim.  These  are  to  be  obtained  through  experimental  photographs. 
The  various  fabric  based  relationships  do  not  account  for  the  micromotions  of  the 
individual  particles  and  thus  is  limited  in  it's  capability  to  predict  the  changes  in 
packing  structure  under  a  deformation.  This  led  to  the  development  of 
microstructural  continuum  theories  which  incorporate  not  only  the  various  fabric 
measures  but  also  the  microstructural  deformations  and  the  micromotions  of  the 
individual  particles. 

2.7  Microstructural  Continuum  theory 

Microstructural  continuum  theory  (Bathurst  and  Rothensburg,  1988;  Chang 
and  Misra,  1989a,b)  is  an  analytical  approach  and  is  implemented  using  computer 
codes  through  incremental  loading.  In  this  theory,  the  micromotions  of  the 


41 


particles  are  expressed  as  linear  or  higher  order  functions  of  the  coordinates  of 
the  centroid  of  particles  and  the  macroscopic  strain  simulated  according  to  the 
type  of  deformation  (biaxial,  simple  shear,  etc.)-  The  theory  incorporates  the 
terms  related  to  the  microstructural  properties  and  also  the  terms  related  to 
packing  structure  such  as  contact  normal,  branch  vector  etc.  The  microstructural 
continuum  theory  was  pursued  as  the  computational  time  required  was  less  than 
that  taken  by  the  computer  simulation  methods  such  as  DEM  and  DMA.  Both  the 
microstructural  continuum  theory  and  the  DEM  idealize  a  granular  assembly  as 
an  assembly  of  rigid  disks/ spheres  linked  by  normal  and  tangential  springs  at 
contacts  (Fig.  2.2).  The  difference  between  the  two  approaches  is  in  the 
transmission  of  boimdary  disturbances  to  the  internal  particles.  In  the  DEM,  the 
boundary  displacements/ forces  are  transferred  to  the  neighboring  particles  one 
by  one  and  the  particles  are  brought  to  equilibrium  through  an  iterative  process. 
In  the  micro-structural  continuum  method,  the  particle  displacements  within  a 
particulate  assembly  are  assumed  to  be  a  polynomial  function  of  the  macroscopic 
strain  applied  at  the  boimdary  and  the  particle  coordinates.  Thus,  the 
macroscopic  strains  are  related  to  particle  displacements  which  in  turn  are  related 
to  the  displacements  of  the  contacts.  The  contact  displacements  are  obtained  from 
the  components  of  the  relative  displacements  of  the  particles,  as  the  particles  are 
assumed  to  be  rigid.  The  contact  displacements  are  related  to  the  contact  forces 
through  the  Hertz-Mindlin  contact  theory  which  in  turn  are  related  to  the 
macroscopic  stress  through  Hill's  (1963)  averaging  principle  described  in  Chapter. 
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4.  Thus,  a  constitutive  relationship  is  established  between  the  macroscopic  stress 
and  macroscopic  strain  as  a  function  of  the  microstructural  properties  and 
packing  structure.  The  packing  structure  is  also  updated  at  the  end  of  an 
incremental  deformation  according  to  the  particle  locations  and  the  loss  and  gain 
of  contacts. 

Bathurst  and  Rothensburg  (1988)  derived  explicit  expressions  for  elastic 
constants  in  terms  of  microstructure  and  verified  the  expressions  with  the  results 
from  numerical  tests.  In  their  study,  the  micromotions  of  particles  were 
incorporated.  However,  the  contacts  were  assumed  to  be  indestructible  and  thus 
the  evolution  of  packing  structure  with  deformation  was  not  considered. 

Chang  and  Misra  (1989a)  evolved  stress-strain  relationships  for  regular 
cubical  and  hexagonal  packings  of  disks  based  on  microstructural  continuum 
theory  and  obtained  good  comparison  with  laboratory  biaxial  tests  at  very  low 
strains  of  =  0.00015.  Chang  and  Misra  (1989b)  developed  a  constitutive 
relationship  assuming  a  linear  displacement  field  for  the  particles.  In  this  model, 
the  evolution  of  packing  structure  for  a  random  assembly  of  particles  during  a 
deformation  was  also  included.  The  results  from  the  model  had  good  agreement 
with  the  numerical  tests  conducted  using  the  numerical  model  of  Serrano  and 
Rodriguez-Ortiz  (1973)  at  low  strain  levels  of  =  0.0005  in  biaxial  tests  on 
assemblies  of  rods.  Constitutive  relationships  incorporating  second  order 
polynomial  functions  for  particle  displacements  and  particle  rotations  can  be 
foimd  in  Chang  and  Liao  (1990),  Chang  (1990)  Chang  and  Misra  (1990)  and  Misra 


43 


(1990). 

An  experimental  effort  to  provide  data  to  relate  the  macroscopic  strain  to 
the  micromotions  of  the  particles  was  carried  out  by  Gill  (1993)  on  granular 
materials  under  uniaxial  strain.  The  conclusion  was  that  the  displacement  fields 
are  not  directly  proportional  to  applied  macroscopic  strain  as  predicted  by  local 
theories.  Strain  in  the  immediate  vicinity  of  the  loading  cap  was  foimd  to  be 
greater  than  the  applied  macroscopic  strain  with  other  regions  of  the  specimen 
generally  experiencing  strain  levels  less  than  the  applied  macroscopic  strain.  Gill 
(1993)  observed  that  displacements  within  a  heighborhood  of  grains  are  non- 
uniform  and  obeys  an  exponential  distribution  and  proposed  a  theory  based  on 
the  diffusion  equation.  However  no  theoretical  relationships  were  derived  for  the 
micromotions  of  particles  based  on  the  diffusion  theory.  The  experimental  study 
established  the  fact  that  displacement  fields  are  not  directly  proportional  to  the 
macroscopic  strain,  a  fundamental  assumption  in  the  microstructural  continuum 
theories. 

2.8  Lattice  model 

A  lattice  approach  coupled  with  a  statistical  analysis  has  been  used  by 
physicists  to  study  the  failure  of  brittle  materials  (Krajdnovic  and  Silva,  1982; 
Hermann  and  Roux,  1990).  Historically,  first  application  of  lattices  to  structural 
problems  is  credited  to  Hrenikoff  (1941).  In  the  lattice  approach  in  which  a 
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material  is  simulated  by  a  system  of  bars /links,  was  imdertaken  from  the  fact  that 
a  lattice  represents  a  physically  justified  and  computationally  efficient 
approximation  of  a  continuum.  Great  advantages  of  the  lattices  was  that  the 
analyses  of  dilution  processes  (sequential  removal  of  links)  and  the  determination 
of  the  influence  of  the  disorder  (introduced  by  a  random  spatial  distribution  of 
the  link  parameters)  are  computationally  not  very  demanding  (Hermann  and 
Roux,  1990).  The  important  features  of  lattices  are:  (a)  ability  to  consider  more 
complex  states  of  stresses  and  strains  and  (b)  efficiency  in  providing  estimates  of 
the  influence  of  stress  concentrations  and  spatial  correlations  of  defects  on  the 
macro-response  of  the  system. 

Krajcinovic  and  Silva  (1982)  used  a  parallel  bar  model  to  study  the  failure 
of  a  perfectly  brittle  material  subjected  to  tension.  In  this  approach,  the  material 
was  replaced  by  a  -system  of  parallel  bars  of  identical  stiffnesses  connected  to 
hinges  at  one  end  and  to  a  very  rigid  member  at  the  other  end.  All  bars  were 
assumed  to  be  perfectly  elastic  and  an  initial  disorder  was  introduced  by 
assuming  the  rupture  strengths  of  the  individual  links  to  follow  a  probability 
distribution  function.  The  elongation  of  the  system  was  then  expressed  as  a 
function  of  the  applied  force  and  a  parameter  called  "damage",  where  damage 
was  defined  as  a  ratio  between  the  number  of  ruptured  bars  and  the  total  number 
of  bars.  The  extent  of  damage  in  turn  was  a  function  of  the  rupture  strength 
distribution  and  the  applied  force.  This  model  was  useful  in  predicting  the 
behavior  of  certain  classes  of  brittle  ceramics  reinforced  with  continuous  fibers 
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and  long  cortical  bones.  A  limitation  of  this  model  was  that  it  ignored  the  stress 
fluctuations  and  defect  interaction  and  this  led  to  the  development  of  models 
consisting  of  triangular  lattices  (Hansen  et.  al  1989).  Thus,  the  emphasis  of 
physicists  in  using  lattices  was  for  the  study  of  disorder  in  materials  through 
statistical  methods  (Krajcinovic  and  Basista,  1991). 

The  interests  of  engineers  have  been  to  evolve  continuum  like  models  for 
materials  with  lattice-like  microstructure  or  to  use  lattices  for  the  analysis  of 
continuum.  The  relationship  between  a  continuum  and  a  gridwork  of  discrete 
elements  was  first  examined  in  some  early  work  by  Hrennikoff  (1941)  and  later 
by  Newmark  (1949)  which  involved  analyzing  a  continuum  by  replacing  it  with 
an  equivalent  elastic  gridwork.  The  reverse  possibility  of  replacing  a  large 
repetitive  gridwork  with  an  equivalent  continuum  was  studied  by  Kollar  et.  al 
(1985),  Dow  et.al  (1985)  and  Noor  and  Russell  (1986).  These  studies  presented 
methods  to  simplify  certain  calculations  for  large  space  trusses  and  space  frames. 
Recently  Dai  and  Frantziskonis  (1994)  used  the  lattice  approach  to  model  a  brittle 
material  through  the  data  obtained  from  ultrasonic  scanning.  The  data  was  used 
to  assign  failure  stress  and  stiffness  for  each  bar  in  the  discretized  members  of  the 
lattice.  The  lattice  model  predicted  crack  patterns  similar  to  that  were  observed 
in  experiments.  Also,  Frantziskonis  et.al  (1994)  used  a  lattice  approach  to  study 
the  fiber  matrix  interface  in  composite  materials. 

A  lattice  approach  called  discrete  stiffness  model  (DSM)  was  developed  by 
Trollope  and  Burman  (1980)  to  study  the  mechanics  of  granular  media.  The 
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model  was  originally  developed  for  the  movement  of  blocks  of  rocks  and  was 
later  extended  for  granular  media.  The  lattice  arrangement  was  formed  by 
linking  the  particle  centroids  as  nodes  with  bars  having  a  joint  element  between 
the  nodes  and  at  the  contact  point.  The  details  of  the  lattice  model  and  joint 
element  are  elaborately  dealt  in  Burman  (1974)  and  are  not  reviewed  here  as  this 
would  involve  extensive  discussion  for  a  proper  imderstanding  of  the  approach. 

The  discrete  stiffness  model  was  applied  on  granular  wedges  to  assess  the 
stress  distribution  and  stability  of  the  slopes  under  gravity  load.  The  study  was 
chosen  on  granular  wedges  as  they  have  direct  relevance  to  practical  structures 
such  as  earth  and  rock-fill  dams,  road  and  rail  embankments,  etc.  The  normal 
stress  distribution  on  the  bottom  boimdary  from  the  model  had  good  agreement 
with  experimental  observation.  The  DSM  was  used  to  study  the  variation  of 
normal  and  shear  stresses  at  the  base  of  granular  wedges  and  the  displacement 
patterns  due  to  vertical  settlements  of  the  wedges  simulated  by  downward 
displacements.  The  main  conclusion  of  the  study  was  that  the  displacement 
patterns  do  not  reflect  the  internal  stress  distribution.  The  study  suggested  that 
interpretation  of  displacement  monitoring  of  earth  and  rock  structure  requires 
careful  evaluation.  In  the  DSM,  the  contacts  were  considered  to  be  indestructible 
and  thus  evolution  of  the  packing  structure  was  not  considered  and  the  study 
was  confined  to  granular  wedges  imder  gravity  load. 

Following  an  approach  similar  to  DSM,  Boardman  (1990)  and  Sadd  et.  al 
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(1992)  studied  the  wave  propagation  in  granular  media.  The  lattice  was  formed 
by  linking  particle  centroids  and  elements  were  modeled  as  1-D  d5mamic  bar 
element  with  lumped  mass  at  the  end  nodal  points  (Fig.  2.7a).  The  element 
possessed  both  normal  and  tangential  stiffness  and  damping  coefficients  whose 
numerical  values  were  obtained  from  dynamic  photoelastic  experiments.  A 
modified  bar  element  incorporating  the  contact  characteristics  was  used  (Fig.  2.7). 
In  this  method,  there  was  no  joint  element  as  in  DSM.  The  degree  of  freedom 
at  a  node  consisted  of  the  normal  and  tangential  displacements  and  a  rotation 
(rotation  of  particles).  The  stiffness  matrix  was  constituted  from  the  normal 
stiffness  {K^),  tangential  stiffness  (K^)  and  the  radii  of  the  particles  (Figs.  2.7b 
and  2.7c). 

The  finite  element  model  was  applied  on  ID  and  2D  assemblies  of  disks 
made  of  Homalite-100  subjected  to  a  dynamic  load  on  a  disk  at  the  top  of  the 
assemblies.  The  resulting  wave  speed  and  the  magnitude  and  orientation  of  the 
contact  forces  obtained  from  the  lattice  model  had  good  agreement  (10-20% 
accurate)  with  the  results  from  photoelastic  experiments.  However,  the 
application  of  the  lattice  model  was  limited  to  small  deformations  and  thus  loss 
and  gain  of  contacts  due  to  fabric  changes  were  not  considered. 


Disk  1 


(a)  Basic  element  of  the  elastic  network 


Kn  -  normal  stiffness 
Kt  —  tangential  stiffness 
JJ.Ta  -  radii  of  particles 


(b)  Disk  equilibrium  for  (c)  Disk  equilibrium  for 

normal  forces  tangential  forces 


Fig.  2.7  Elastic  network  modeling  of  granular  media  (Sadd  et.al,  1992) 
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In  particulate  materials,  loads  are  transmitted  along  discrete  paths  within 
the  material  through  particle  contacts.  The  dispersion  of  loads  witiiin  the  material 
can  be  envisioned  to  be  through  a  network  formed  by  linking  the  particle 
contacts.  It  is  this  network  pattern  of  load  distribution  within  a  particulate 
material  that  gave  rise  to  the  idea  of  lattice  type  simulation  of  particulate 
materials  presented  in  this  research.  The  lattice  simulation  is  carried  out  by 
linking  the  particle  contacts  as  different  from  the  lattice  approaches  of  Trollope 
and  Burman  (1980)  and  Sadd  et.  al  (1992)  in  which  the  lattice  was  formed  by 
linking  the  particle  centroids.  Thus  the  lattice  type  model  in  this  research 
attempts  to  model  a  particle  itself  as  a  lattice  within  the  microstructure  of  the 
particle. 


CHAPTER  3 


FORMULATION  OF  THE  LATTICE  TYPE  MODEL 
FOR  PARTICULATE  MEDIA 


3.1  Model  Formulation 


3.1.1  General 


In  the  lattice  type  simulation  of  a  two-dimensional  assembly  of  particles, 
the  assembly  is  transformed  into  a  network  of  bars  representing  the  load  path 
within  the  particles  (Fig.  3.1).  Nodes  are  introduced  at  particle  contacts  and  are 


Particle  (typ) 


Bar  (typ) 


I  kvWW 
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Node  (typ) 


a)  (b) 

Fig.  3.1  (a)  Particle  arrangement  (b)  Lattice  type  simulation 
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linked  by  bars.  A  contact  node  on  a  particle  is  linked  to  every  other  contact  node 
on  that  particle  with  bars  resulting  in  a  truss.  If  the  number  of  contacts  around 
a  particle  is  less  than  four,  then  bondless  nodes  are  introduced  around  the 
particle,  so  as  to  have  a  minimum  of  four  nodes  forming  the  lattice  for  the 
particle.  This  is  necessary  because,  during  the  deformation  process,  the  trajectory 
of  the  centroid  of  a  particle  is  traced  from  the  displacements  of  the  centroid  of  the 
polygon  formed  by  linking  the  contacts  around  that  particle.  Particles  of  irregular 
shape  can  be  approximated  to  disks  of  circular/elliptical  cross-sections  or  mapped 
using  the  theory  of  morphology  (Luerkens,  1991).  However,  this  is  not  pursued 
at  this  stage  of  the  research. 

The  elastic  modulus  of  the  bars  within  a  particle  is  assumed  to  be  equal  to 
the  elastic  modulus  of  the  particle.  The  lengths  of  the  bars  depend  on  the 
geometry  of  the  particles  and  the  location  of  the  contact  points.  Since  the  bar 
areas  are  dependent  on  the  bar  forces  which  are  not  known  initially,  the  initial 
cross-sectional  areas  of  all  the  bars  in  the  lattice  are  proportioned  to  the  initially 
applied  boundary  load  (P)  on  a  boundary  particle. 

Consider  the  circular  disk  imder  a  diametrically  opposite  load  in  Fig.  3.2. 
The  initial  cross-sectional  area  (A)  of  the  bars  in  the  lattice  is  set  equal  to  the  area 
of  a  bar  of  length  equal  to  the  average  diameter  of  the  particles  (L)  and  subjected 
to  a  load  P  that  will  give  the  displacement  (6)  in  the  above  elastic  disk.  'A'  is 
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given  by  A  =  (PL)/ (E6);  where  E  is  the  elastic  modulus  of  the  particle  and  6  is 
the  relative  vertical  displacement  of  A  relative  to  C  in  the  disk  shown  in  Fig.  3.2. 
With  these  parameters  (lengths,  areas  and  elastic  moduli),  the  truss  is  analyzed 
for  a  given  load 

.  P 


P 

Fig.  3.2  Disk  imder  a  Vertical  Load 

imder  appropriate  boundary  conditions  using  the  "standard"  governing 
equation  (3.1)  in  structural  mechanics: 

[K]  {0}  =  {R)  (3.1) 

where  [K]  is  the  global  stiffness  matrix  of  the  truss,  {D}  is  the  displacement  vector 
comprised  of  the  displacements  of  the  nodes  of  the  truss  in  the  X  and  Y 
coordinate  directions  and  {R}  is  the  global  load  vector  comprised  of  the  loads 
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acting  at  the  nodes  of  the  truss  in  X  and  Y  directions.  [K]  is  a  square  matrix  of 
size  "2Nx2N"  where  'N'  is  the  total  number  of  nodes  in  the  truss  structure  and 
{D}  and  {R}  are  column  vectors  of  size  "2Nxl".  The  above  equation  is  solved  for 
the  imknown  displacement  vector. 

3.1.2  Formulation  of  Global  Stiffness  matrix 

The  lattice  t57pe  model  (LTM)  is  a  global  technique  which  simulates  the 
entire  assembly  of  particles  in  a  particulate  medium  as  a  truss  structure. 
Therefore,  the  first  step  in  the  computational  process  is  to  evaluate  the  global 
stiffness  matrix  [K].  The  global  structure  stiffness  matrix  [K]  is  generated  by 
taking  the  sum  of  element  stiffness  matrices  of  the  truss  (bar)  elements  using 

nel 

[K]  =  E  W  (3.2) 

/J=:1 

where  [k]  is  the  element  stiffness  matrix  and  "nel"  is  the  total  number  of  elements 
in  the  truss  structure.  The  element  stiffness  matrix  [k]  for  a  truss  element  is  a 
function  of  the  elastic  modulus,  length  and  area  of  cross-section  of  the  bar 
element  and  also  the  sines  and  cosines  of  ti\e  angle  made  by  the  orientation  of  the 
bar  with  the  positive  direction  of  X-axis  (Fig.  3.3a). 

Consider  a  bar  element  of  length  'U  and  cross-sectional  area  'A'  shown  in 


5' 


Fig.  3.3a  with  a  back  node  'i'  and  fore  node  'j'  and  which  makes  an  angle  p  with 
the  positive  direction  of  X-axis.  Let  the  elastic  modulus  of  the  material  of  the  bar 
be  'E'.  The  input  data  for  the  analysis  are  the  X  and  Y  coordinates  of  the  nodes 
of  the  bar  'i'  and  'j'.  The  matrix  [k]  is  generated  by  computing  the  length  and  the 
sines  and  cosines  of  P  from  Equations  (3.3)  and  (3.4)  respectively. 


(a)  arbitrarily  oriented  (b)  after  imposed  displacement 

in  the  xy  plane  v,  =  Uj  =  Vj  =  0 


Fig.  3.3  Truss  element 

L  =  [(Xj  -  xf  +  {yj  -  (3.3) 


c 


COSp  = 


(3.4) 
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The  degrees  of  freedom  in  a  truss  are  the  X  and  Y  displacements  of  the 
nodes.  The  displacements  along  the  degrees  of  freedom  for  the  bar  in  Fig.  3.3  are 
denoted  by  Uj,  Vj  at  node  'i'  and  Uj,  Vj  at  node  'j'.  The  elements  of  tibe  matrix  [k] 
are  generated  by  assembling  the  column  vectors  (matrices)  obtained  by  imposing 
a  displacement  in  the  direction  of  each  degree  of  freedom  (d.o.f)  while  other 
degrees  of  freedon  are  set  at  zero.  The  first  of  these  four  cases  is  shown  in  Fig. 
3.3b  in  which  a  displacement  is  imposed  at  node  'i'  in  the  X-direction.  This 
displacement  induces  an  axial  shortening  equal  to  CUj  of  the  bar  element  which 
in  turn  produces  an  axial  compressive  force  'F'  given  by 

F  =  (AE/L)ci7,.  (3.5) 

The  X  and  y  components  of  F  at  nodes  'i'  and  'j'  are  given  by  Pj  =  -pj  =  Fc  and 
Qj  =  -Pj  =  Fs.  These  force  components  Pj,  Pj,  Pj  are  grouped  as  a  column 
vector  and  is  called  the  load  vector  for  the  bar  element.  When  the  expression  for 
'F'  in  Equation  (3.5)  is  substituted  in  the  expressions  for  P/,  Pj,  Pj,  Pj,  there  is  a 
common  factor  {AE/L)  u,..  Taking  this  common  factor  out  of  the  expressions  for 
P/I  Pj,  Pj,  Pj,  the  load  vector  is  expressed  in  a  matricial  form  as 
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In  a  similar  manner,  displacements  can  be  imposed  along  the  other  degrees  of 
freedom  as  Vj,  Uj  and  Vj  and  equations  similar  to  Equation  (3.6)  can  be  obtained 
for  each  of  these  degrees  of  freedom.  Superimposing  the  four  equations  for  the 
four  degrees  of  freedom.  Equation  (3.7)  is  obtained,  which  represents  the  element 
level  equation  for  a  bar  subjected  to  imposed  displacements  of  Uj,  Vj  at  node  'i' 
and  Uj,  Vj  at  node  'j'. 
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(3.7) 


The  term  AE/L  multiplied  by  the  square  matrix  in  Equation  (3.7)  is  referred  to 
as  the  element  stiffness  matrix  [k].  Representing  the  vector  comprised  of  the 
displacement  componets  as  the  displacement  vector  {d}  and  the  load  vector  at 
element  level  as  {r}.  Equation  (3.7)  is  rewritten  as  Equation  (3.8)  which  is  called 
the  equilibrium  equation  of  the  element. 

m  {d}  =  {!)  (3.8) 

Summation  of  the  element  stiffness  matrices  [k]  results  in  the  global  stiffness 
matrix  [K]  given  by  Equation  (3.2). 
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3.1.3  Formulation  of  Global  load  vector 

There  are  three  t5^es  of  loads  imposed  on  a  truss  structure  and  each  of 
these  is  incorporated  in  the  global  load  vector  {R}  by  different  computational 
procedures.  The  three  types  of  loads  are 

(1)  Through  prescribed  nodal  loads 

(2)  Through  prescribed  nodal  displacements 

(3)  Through  prescribed  internal  initial  axial  forces  in  the  bars. 

In  the  previous  section,  it  was  shown  through  Equation  (3.7)  that  if  the 
displacements  are  prescribed,  then  the  resulting  forces  at  the  ends  of  the  bar  to 
keep  it  in  equilibrium  can  be  determined.  Contrary  to  this,  if  the  loads  are 
prescribed,  then  the  displacements  can  be  determined.  Unlike  the  case  of  a  single 
bar  element,  in  a  global  truss  structure,  the  displacement  or  a  load  at  a  node  will 
have  global  influence.  The  computational  procedures  used  to  include  the  effects 
of  the  above  loading  t5q?es  are  explained  here.  The  elements  in  the  load  vector 
{R}  corresponding  to  node  'j'  are  denoted  as  R  {2j-^)  and  R  (2j),  which 
represents  the  loads  in  the  X  and  Y  directions  respectively.  The  elements  of  the 
load  vector  {R}  are  initially  set  to  zero  and  then  assembled  to  include  the  effects 
of  various  types  of  loads. 

Let  andPj  be  the  prescribed  nodal  loads  in  a  truss  at  a  node  'j'  in  the 
X  and  Y  directions.  Then  the  effects  of  these  are  included  in  the  load  vector  {R} 
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by  adding  the  loads  to  the  components  of  the  load  vector  algebraically  as 
R  (2y-1)  =  R  (2y-1)  +  and  R  (2y)  =  R  {2j)  +  P^. 

The  second  loading  condition  of  prescribed  nodal  displacements  is 
incorporated  in  the  load  vector  by  first  modifying  the  load  vector  and  then  the 
stiffness  matrix.  The  reader  may  refer  to  Cook  (1981)  for  details  additional  to  the 
ones  presented  here.  Consider  a  stipulated  displacement  6  at  a  node  'k'  in  a  truss 
in  the  X-direction,  whose  d.o.f  number  is  denoted  as  n  =  (2k  -  1).  The  effect  of 
this  is  included  in  the  load  vector  by  modifying  the  load  vector  {R}  through 
Equation  (3.9)  for  every  element  R(i)  of  the  load  vector,  with  'i'  ranging  from  1 
to  2N,  where  N  is  the  total  number  of  nodes  in  the  truss. 

R{i)  =  R{i)  -  Km  6  (3.9) 

After  modifying  the  load  vector  for  all  the  prescribed  nodal  displacements,  the 
load  vector  is  again  modified  by  setting  the  elements  of  the  load  vector  equal  to 
the  prescribed  displacements  only  for  d.o.f  where  displacements  are  prescribed. 
Therefore,  for  the  above  example,  the  load  vector  is  modified  as  R(n)  =  6.  After 
the  above  two  procedures,  the  stiffness  matrix  is  modified  by  setting  K(n,n)  =  1 
and  all  the  other  elements  of  the  row  'n'  and  column  'n'  of  [K]  equal  to  zero. 

The  effect  of  initial  internal  forces  in  the  bars  of  a  truss  are  incorporated 
in  the  load  vector  {R}  by  adding  the  components  of  the  forces  at  the  two  ends  of 
the  bars  to  the  components  of  {R}  corresponding  to  the  nodes.  An  internal  force 
'F  in  a  bar  shown  in  Fig.  3.3  exerts  a  force  of  Pj  =  -pj  =  Fc  and  q.  =  -Qj  =  Fs 
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which  is  expressed  in  a  matricial  form  as  Equation  (3.10). 
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Then,  the  elements  of  the  load  vector  {R}  is  modified  as  follows. 
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3.1.4  Computation  of  Axial  forces 

After  evolving  the  global  stiffness  matrix  and  the  global  load  vector,  the 
next  step  is  the  solution  of  Equation  (3.1)  which  results  in  the  displacement  vector 
{D}.  The  solution  is  carried  out  using  the  Gaussian  elimination  method.  From 
the  displacement  of  the  nodes,  the  components  of  the  load  acting  at  the  two  ends 
of  an  element  are  computed  using  Equation  (3.7).  From  these  components,  the 

Pi 

axial  force  in  a  bar  'k'  is  computed  as  =  —  . 
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3.1.5  Normal  force.  Shear  force  and  Shear  force  ratio 

When  two  particles  are  pressed  together  with  a  contact  force,  a  contact 
surface  develops  and  the  contact  normal  is  defined  as  the  direction  normal  to  the 
contact  surface  (Fig.  2.5b).  For  two  particles  of  circular  cross-sections,  the  contact 
normal  is  the  direction  linking  the  centroids  of  the  two  particles  in  contact  (Fig. 
3.4). 

The  normal  force  and  shear  force  at  a  contact  node  are  computed  by  taking 
the  algebraic  sum  of  the  components  of  the  axial  forces  in  the  bars  in  the  direction 
of  the  contact  normal  and  tangential  to  the  contact  normal  (Fig.  3.4)  and  are  given 
by 


Contact  normal 
Shear  force 
Normal  force 

Resultant  forces 
at  contact  nodes 


Fig.  3.4  Particle  under  contact  forces 
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n~1 

=  E  sin6y  (3.12) 

y=i 

S,  .  S  Fj  cos9^  (3.13) 

y=i 

where  Nj  =  normal  force  at  the  contact  node  'i';  Sy  =  shear  force  at  the  contact 
node  'i';  Fj  =  axial  force  in  bar  6j  =  angle  made  by  bar  with  the  shear 
force  direction  at  the  contact  node  'i';  n  =  number  of  contacts  around  a 
particle.  Shear  force  ratio  at  a  contact  node  T  denoted  by  A}  is  the  ratio  of  the 
absolute  value  of  the  shear  force  divided  by  the  normal  force  and  is  given  by 

=  J  (3.14) 

/ 

Since  the  contact  area  for  the  computation  of  the  contact  shear  stress 
and  contact  normal  stress  are  the  same,  the  shear  stress  ratio  defined  as  the 
shear  stress  divided  by  the  normal  stress  is  equal  to  the  shear  force  ratio. 
Therefore,  the  shear  force  ratios  will  also  be  referred  to  as  shear  stress  ratios. 

3.1.6  Cross-sectional  areas  of  bars 

The  cross-sectional  areas  of  the  bars  are  updated  at  the  end  of  each  load 
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increment  based  on  the  new  contact  forces  around  the  particles.  The  area  of  a  bar 
is  found  by  computing  the  displacements  of  an  elastic  disk  under  contact  forces 
and  then  determining  the  bar  area  that  will  give  an  equivalent  displacement  in 
the  direction  of  the  bar.  In  this  research,  only  particles  with  circular  cross-sections 
are  considered  and  thus  only  the  method  to  compute  the  displacements  in  a 
circular  disk  is  explained  here.  However,  the  method  can  be  extended  to  particles 
of  elliptical  or  other  shaped  cross-sections. 

The  state  of  stress  in  the  two  dimensional  element  of  a  continuum  is  shown 
in  Fig.  3.5.  The  stresses  o^,  0^  and  are  the  normal  stress  in  the  X-direction, 
normal  stress  in  the  Y-direction  and  shear  stress.  From  theory  of  elasticity. 


Fig.  3.5  State  of  stress  in  the  two  dimensional  element 
of  a  continuum 
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Michell  (1900)  derived  expressions  for  the  stresses  o^,  Oy  and  T^^y  in  an  infinite 
half  space  of  continuum  acted  upon  by  a  point  load.  From  these  expressions, 
Michell  (1902)  derived  expressions  for  O^,  Oy  and  T^y  at  a  point  in  a  circular  disk 
acted  upon  by  an  equal  and  opposite  load  acting  along  a  chord  of  the  disk.  The 
displacements  of  the  points  of  a  circular  disk  acted  upon  by  a  set  of  contact 
forces  is  derived  from  the  above  expressions  for  stresses.  The  circular  disk 
problem  and  the  derivation  of  the  stresses  are  also  given  in  Timoshenko  and 
Goodier  (1970)  and  in  Frocht  (1948). 

The  X  and  Y  components  of  the  displacements  of  the  contact  points  of  the 
disk  shown  in  Fig.  3.4  are  determined  by  superimposing  the  displacements 
produced  by  each  contact  force  (Fig.  3.6).  Consider  only  the  contact  force  PI 
acting  on  the  disk  (Fig.  3.6b)  at  node  1.  Then,  the  resultant  of  P2,  P3  and  P4 
denoted  as  P5  and  equal  to  PI  in  magnitude  will  be  acting  at  node  5  along  the 
line  of  action  of  PI  and  opposite  to  the  direction  of  PI.  The  displacements 
produced  by  PI  and  P5  will  be  symmetrical  about  the  line  AB,  a  diameter  of  the 
circle  normal  to  the  line  of  action  of  PI.  Then,  it  can  be  seen  that  the 
displacements  produced  by  the  load  PI  on  the  circular  disk  is  only  in  the  semi¬ 
circular  portion  of  the  disk  towards  the  right  of  AB.  Therefore,  the  influence  of 
PI  is  only  at  nodes  1  and  4.  The  displacements  due  to  load  PI  are  equal  to  6.,., 
at  node  1  and  at  node  4.  The  displacements  6.,.,  and  64.,  are  obtained  from 
the  expressions  for  displacements  in  a  circular  disk  acted  upon  by  an  equal  and 
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+ 

AB  -  Diameter  normal  to 
line  of  action  of  PI 


d..  =  displacement  of  node  i 
due  to  load  at  node  j 


Aj  =  net  displacement  at  node  j 


Fig.  3.6  Superposition  of  Displacements 


opposite  load  acting  along  a  chord  of  the  circle  similar  to  the  disk  in  Fig.  3.6b. 
The  expressions  for  the  displacements  are  derived  in  Appendix  -  1. 

Following  the  above  approach,  the  displacements  caused  by  the  loads  P2, 
P3  and  P4  acting  at  nodes  2,  3  and  4  at  their  points  of  application  and  at  other 
nodes  are  also  determined.  Finally,  the  resultant  displacements  A.,,  Ztj,  Ag  and  A^ 
at  nodes  1,  2,  3  and  4  are  determined  using  the  rule  of  superimposition  as 


Ai 


=  6 


11 


(3.15) 
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Aj  -  522 

(3.16) 

^3  “  ®33 

+  ^32 

(3.17) 

A4  =  644 

5^3 

+  5..  +  — 

41  2 

(3.18) 

The  displacements  of  the  contacts  are  resolved  in  the  direction  of  each  bar 
and  the  relative  displacements  in  the  direction  of  each  bar  are  computed.  The  bar 
areas  are  then  updated  using 


E,  6,. 


(3.19) 


where  Aj  =  cross-sectional  area  of  bar  'j';  Fj  =  axial  force  in  the  bar;  Lj  =  length 
of  the  bar;  Ej  =  Elastic  modulus  of  the  bar  and  5y  =  relative  displacement  of  the 
nodes  of  bar  in  the  direction  of  bar  'f. 

As  an  example  of  a  lattice  type  simulation,  consider  a  circular  disk  (Fig. 
3.7)  under  a  diametrically  opposite  load  and  restrained  by  lateral  supports.  The 
vertical  displacement  at  A  can  be  found  using  the  expressions  for  displacements 
in  Appendix  -  1.  Let  6^^  be  the  horizontal  displacement  at  D  due  to  the  vertical 
load,  if  there  were  no  lateral  supports.  If  a  horizontal  load  'R'  acts  at  B  and  D  in 
opposite  direction  which  will  produce  a  horizontal  displacement  (5^^)  equal  to 
,  then  'R'  will  be  the  lateral  reaction  in  the  disk  restrained  by  lateral  supports. 


Fig.  3.7  (a)  Disk  under  vertical  load  (b)  Equivalent  lattice 


The  expression  for  5^^  is  given  in  Appendix  - 1,  which  will  be  a  function  of  the 
lateral  support  reaction  'R'.  By  setting  equal  to  6^^,  'R'  can  be  found.  By 
considering  the  force  equilibrium  equations  at  A  and  D,  the  bar  forces  can  be 
evaluated.  Resolving  the  displacement  of  A  in  the  direction  of  bars  AC  and  AD, 
the  deformations  of  these  bars  can  be  evaluated.  Then  the  equivalent  bar  areas 
can  be  computed  from  Equation  (3.19). 

If  a  bar  is  in  tension,  the  axial  force  in  that  bar  is  set  to  zero  and  the 
unbalanced  forces  are  redistributed  through  an  initial  stress  analysis.  The  area 
of  the  bar  is  then  reduced  to  a  low  value  equal  to  10'®  times  the  maximum  of  the 
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cross-sectional  area  of  all  the  bars  within  the  lattice  forming  the  particulate 
assembly.  The  acronym  "AWBAR"  is  used  to  denote  this  reduced  area. 

Since  the  cross-sectional  areas  of  all  the  bars  are  set  equal  to  the  same  value 
for  the  initial  load  analysis,  equivalent  bar  areas  which  will  simulate  the  exact 
contact  force  distribution  within  a  particulate  assembly  under  the  initial  load  are 
accomplished  through  an  iterative  process.  The  bar  areas  are  updated  after  the 
analysis  imder  the  initial  load.  Again  the  initial  load  analysis  is  carried  out  with 
the  updated  bar  areas  and  again  the  bar  areas  are  updated.  The  process  is 
repeated  imtil  the  maximum  bar  area  within  the  lattice  converges  to  an  accuracy 
of  six  decimal  places  with  respect  to  the  bar  area  in  previous  iteration  step.  At 
this  stage,  additional  load  is  applied  on  the  particulate  assembly. 

3.1.7  Redistribution  of  Additional  shear  forces 

If  the  shear  stress  ratio  at  a  contact  node  exceeds  the  coefficient  of  friction 
of  the  material  (p),  the  additional  shear  force  is  redistributed  to  the  other  contacts 
by  applying  an  equal  and  opposite  shear  force  at  the  contact  as  shown  in  Fig.  3.8 
and  given  by  Equation  (3.20). 

AS,.  =  (0  -  M)  W,  (3.20) 

where  AS,-  is  the  additional  shear  force  at  the  contact  node  'i',  r,-  is  the  shear 
stress  ratio  and  Nj  is  the  normal  force. 


Hinge 
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Fig.  3.8  Redistribution  of  additional  shear  forces 


Every  bar  has  two  nodes  referred  to  as  fore  node  and  back  node.  When 
an  additional  shear  force  is  redistributed  at  a  contact  node,  then  the  contact  node 
is  designated  as  fore  node  and  the  other  node  as  the  back  node  for  all  the  bars 
linked  to  the  contact  node.  The  additional  shear  force  is  applied  at  the  contact 
node  and  resolved  into  the  bars  as  axial  forces  relative  to  the  stiffnesses  of  the 
bars  and  assuming  the  back  nodes  of  the  bars  as  hinged  (Fig.  3.8).  These  axial 
forces  are  applied  as  initial  axial  loads  in  the  bars  and  an  initial  stress  analysis  is 
carried  out.  Thus,  the  shear  stress  ratios  at  the  contacts  which  exceeded  'p'  are 
reduced.  The  additional  shear  forces  are  redistributed  at  all  the  nodes  where  the 
shear  stress  ratios  exceed  'p'  in  a  single  analysis.  Also,  when  additional  shear 
forces  are  redistributed  at  one  contact  node,  it  may  increase  the  shear  stress  ratios 
at  neighboring  nodes  which  are  above  'p'.  Due  to  the  above  two  reasons,  an 
iterative  process  is  required  to  bring  down  the  shear  stress  ratios  at  the  contact 
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nodes  which  exceeded  to  ']i'. 

3.1.8  Bonding  and  Debonding  of  Contacts 

When  the  normal  force  at  a  contact  becomes  zero,  the  contact  is  split  into 
two  nodes  (debonding).  Constraints  are  imposed  to  check  for  possible  re-bonding 
of  contacts  based  on  the  particle  geometry,  location  of  the  centroids  of  particles 
and  the  equations  (relationship  between  X  and  Y  coordinates)  of  the  imaginary 
plates  representing  the  boimdaries  of  the  particulate  assemblies.  An  example  of 


Fig.  3.9  Bonding  and  Debonding  of  Particles 
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the  bonding  and  debonding  can  be  shown  with  respect  to  the  boundaries  in  a 
simple  shear  apparatus  shown  in  Fig.  3.9. 

Consider  the  particle  'M'  of  radius  'R'  in  the  particulate  assembly  in  Fig. 
3.9  subjected  to  a  shear  strain  'y'.  In  order  to  check  whether  the  particle  'M'  is 
in  contact  with  the  boimdary  AD,  the  following  procedures  are  followed. 

(1)  Derive  the  equation  of  line  AD 

(2)  Draw  a  perpendicular  from  the  centroid  of  particle  to  line  AD 

(3)  Determine  the  perpendicular  distance  GH  as  'D'. 

(4)  If  D  <  (R  +  Ax)  then  the  contact  at  H  is  re-bonded  and  if  D  >  (R  +  Ax)  then 
the  contact  is  de-bonded,  where  Ax  is  a  tolerance  equal  to  0.0001  mm. 

The  centroid  of  the  particle  'M'  is  located  at  (x^  ,  y^).  The  slope  of  the  line 

AD  is,  sp  =  tan  (90  -  tan“\).  Therefore  the  equation  of  AD  is  y  =  sp  x.  Since 


the  product  of  the  slopes  of  two  perpendicular  lines  is  equal  to  -1,  the  slope  of 

<1 

GH  is  - —  from  which  the  equation  of  line  GH  can  be  written  as 
sp 

1 

y  =  - —  X  +  ct,  where  'ct'  is  a  constant.  Since  (x.,  ,  is  a  point  in  GH, 
sp 

1 

/i  =  - —  +  ct  from  which  the  constant  'ct'  can  be  obtained  as 

sp 

/r 

sp 

In  order  to  determine  the  perpendicular  distance  'D',  the  point  of 


intersection  'H'  at  (Xj  ,  /j)  is  to  be  determined.  Since  (Xj  ,  /j)  is  a  point  on  line 
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GH  and  also  on  line  AD,  the  point  (Xj  ,  /j)  can  be  substituted  in  the  equations 

for  lines  AD  and  GH,  derived  above,  and  the  two  equations  can  be  solved  for 

(x.  +  sp  y.) 

Xj,  /j.  This  process  gives  X2  = - and  -  sp  X2.  The  length  of  line 

(1  +  sp2) 

GH  denoted  as  'D'  is  then  determined  as  D  =  ^(x^  -  Xj)^  +  (y^  -  This 
value  of  'D'  is  checked  for  the  criterion  in  step  (4).  Similar  procedures  are 


followed  for  the  right,  top  and  bottom  boundaries. 


When  nodes  split,  the  bar  areas  of  die  bars  linked  to  the  split  nodes  .are  set 
equal  to  AWBAR.  The  nodes  are  renumbered  before  each  analysis  in  order  to 
obtain  a  minimum  bandwidth  for  the  stiffness  matrix. 


3.2  Kinetics  of  particles 

3.2.1  Single  particle  approach 

Particle  kinetics  means  the  sliding,  rotation  and  rolling  of  particles  with 
respect  to  the  neighboring  particles  or  boundaries.  In  the  lattice  type  model,  the 
micro  slips  occurring  at  contacts  are  neglected.  When  the  lattice  equivalent  of  the 
particulate  assembly  in  Fig.  3.1  is  subjected  to  a  load,  the  joints  are  treated  as 
rigid.  Therefore,  during  the  truss  analysis,  relative  movement  between  two 
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particles  due  to  the  sliding  of  particles  at  a  contact  is  not  considered.  A  contact 
node  is  considered  to  have  failed  if  the  shear  force  ratio  at  that  contact  exceeds 
'p'.  A  particle  having  four  contacts  with  only  one  failed  contact  cannot  have 
significant  relative  movement  relative  to  neighboring  particles  or  boundaries. 
However,  if  a  particle  has  'n'  nodes  and  if  all  the  nodes  except  two  nodes  (n-2) 
have  failed,  then  the  particle  can  rotate  and  translate  relative  to  the  neighboring 
particles.  Based  on  this  assumption,  a  method  was  evolved  to  handle  the  kinetics 
of  particles  as  described  below. 

The  rotation  and  translation  of  particles  are  assumed  to  be  caused  by  the 
additional  shear  forces  (Equation  3.20)  acting  at  the  contacts.  The  displacements 
of  the  contact  nodes  due  to  these  forces  are  found  by  carrying  out  an  analysis  of 
the  lattice  within  the  microstructure  of  the  particles  (Fig.  3.10)  by  considering  only 
one  particle  at  a  time.  The  boundary  conditions  in  this  analysis  for  the  bonded 
contacts  are  spring/ elastic  supports.  The  stiffnesses  of  these  springs  are:  normal 
force  divided  by  normal  displacement  and  shear  force  divided  by  tangential 
displacement.  The  loads  are  the  additional  shear  forces  exceeding  the  coefficient 
of  friction  (p)  multiplied  by  the  normal  force.  Tangential  springs  are  not  present 
at  the  nodes  where  additional  shear  forces  are  applied.  The  lattice  within  the 
microstructure  of  the  particles  are  rotated  and  translated  one  by  one,  updating 
the  coordinates  of  the  contacts  each  time.  At  the  end  of  this  process,  an  initial 
stress  analysis  is  carried  out  to  bring  the  truss  in  Fig.  3.1  to  equilibrium  in  the 
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Lattice 


Fig.  3.10  Single  particle  approach  to  particle  kinetics 

updated  coordinate  position. 

After  implementing  the  above  procedure,  the  approach  was  found  to  be 
computationally  inefficient  and  also  erroneous.  The  reasons  for  these  are  to  be 
explained  in  relation  to  some  numerical  values  related  to  a  numerical  test. 
Therefore,  the  simple  shear  test  of  a  two  dimensional  assembly  of  disks 
representing  quartz  grains  shown  in  Fig.  3.11a  is  chosen.  The  simple  shear 
conditions  simulated  here  are  non-ideal  conditions  (Budhu,  1988)  allowing  free 
vertical  movement  of  the  top  boundary.  The  simple  shear  strain  is  applied  on  the 
particle  assembly  by  subjecting  the  lattice  equivalent  of  the  disk  assembly  to  a 
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(b) 


Fig.  3.11  (a)  Development  of  mechanism  in  a  loose  packing  subjected  to  simple 
shear  strain  (b)  Mechanism  framework 
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triangular  variation  of  displacements  as  shown  in  Fig.  3.11a. 

The  particle  assembly  was  subjected  to  incremental  shear  strains  of  10'®. 
This  limiting  value  was  chosen  in  order  to  keep  the  increase  in  the  value  of  the 
shear  stress  ratios  below  20%  of  the  coefficient  of  friction  at  contacts.  With 
incremental  shear  strains  equal  to  10"®,  the  number  of  increments  required  to 
reach  a  shear  strain  of  0.20  is  equal  to  0.20/10"®  =  20,000.  This  is 
computationally  very  inefficient  and  hence  was  found  as  one  limitation  of  the 
single  particle  approach. 

Laboratory  experiments  on  the  loose  packing  similar  to  that  shown  in  Fig. 
3.11a  showed  that  when  the  packing  is  subjected  to  simple  shear  test  conditions, 
the  volume  of  the  packing  decreases  (Dragoo,  1995).  When  the  lattice  equivalent 
of  the  particulate  packing  was  subjected  to  an  incremental  shear  strain,  there  was 
an  increase  in  volume.  When  the  vertical  planes  transformed  into  slippage 
planes,  the  single  particle  approach  was  applied  to  handle  particle  kinetics  and 
this  resulted  in  a  decrease  in  volume.  However,  the  net  change  in  volume  from 
the  analysis  for  incremental  shear  strain  and  from  the  analysis  for  particle  kinetics 
was  an  increase  in  volume  which  is  not  in  accordance  with  the  experiments. 

As  stated  in  Chapter.  1,  predominant  deformation  in  particulate  media  is 
caused  by  the  sliding  and  rotation  of  particles  with  only  a  small  contribution  from 
particle  deformations.  The  above  mentioned  shear  strain  of  10"®  is  in  the  range 
of  the  strains  due  to  deformation  of  particles.  Therefore,  an  alternative  method 
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was  required  to  deal  with  the  kinetics  of  particles.  One  such  alternative  technique 
developed  is  a  global  mechanism  approach  described  in  the  next  section. 

3.2.2  Global  mechanism  approach 

Numerical  and  experimental  observations  on  granular  assemblies  (Cundall 
and  Strack,  1983;  Drescher  and  De  Josselin  De  Jong,  1972)  have  shown  that  large 
forces  are  carried  by  chains  of  particles  and  the  assembly  deforms  like  a  pin- 
jointed  structure  when  a  mechanism  is  reached  and  the  deformation  proceeds 
with  no  increase  in  stress.  Thus,  a  mechanism  in  a  granular  assembly  can  be 
defined  as  a  stage  at  which  the  assembly  cannot  carry  an  additional  load  and  the 
assembly  deforms  under  a  constant  applied  load.  It  was  observed  that  when  a 
granular  assembly  is  subjected  to  a  deviatoric  stress,  particles  lose  contact  in  the 
minor  principal  stress  direction  and  the  number  of  contacts  decreases.  Based  on 
the  above  observations  in  numerical  tests  using  the  DEM,  Cundall  and  Strack 
(1983)  defined  a  term  called  "constraint  ratio"  (CR)  to  check  the  development  of 
mechanism  in  a  granular  assembly  as 

CR  =  (3.21) 


where,  N  =  number  of  particles;  C  =  number  of  contacts;  S  =  fraction  of  sliding 
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contacts  (number  of  sliding  contacts  divided  by  the  number  of  contacts).  The 
numerical  tests  of  Cundall  and  Strack  (1983)  showed  that  when  CR  =  1,  a 
mechanism  had  developed. 

In  the  lattice  type  model,  a  particulate  assembly  is  checked  for  the 
development  of  a  mechanism  after  an  incremental  load  is  applied  on  the  truss 
shown  in  Fig.  3.1.  A  mechanism  is  assumed  to  have  been  reached  if  any  of  the 
following  two  conditions  are  satisfied  (1)  The  increment  in  load  is  zero  due  to 
the  application  of  an  incremental  strain  (such  as  biaxial  or  simple  shear)  on  the 
truss  shown  in  Fig.  3.1.  (2)  The  shear  stress  ratio  at  (n-2)  contacts  around  all  the 
particles  have  reached  'p',  where  'n'  is  the  number  of  contacts  around  a  particle 
while  the  constraint  ratio  is  less  than  1.2.  The  first  condition  is  the  general 
condition  for  a  random  assembly  of  particles.  However,  for  a  regular  assembly 
disks  shown  in  Fig.  3.11a,  the  first  condition  is  not  fulfilled  and  then  the  second 
condition  is  to  be  checked  for  the  development  of  a  mechanism. 

If  a  particulate  assembly  develops  into  a  mechanism,  further  deformation 
is  simulated  through  a  "mechanism  framework"  (Fig.  3.11b).  To  formulate  this 
structure,  nodes  are  introduced  at  particle  centroids  and  are  linked  with  bars. 
Stiffness  of  bars  linking  particle  centroids  with  non-sliding  contacts  are  very  large 
or  "infinite"  (10®)  relative  to  the  stiffness  of  bars  linking  particle  centroids  with 
sliding  contacts.  The  mechanism  analogy  is  explained  with  respect  to  the  simple 
shearing  (non-ideal)  of  an  assembly  of  disks  representing  quartz  grains  shown  in 
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Fig.  3.11a.  In  this  assembly,  a  mechanism  develops  when  all  the  vertical  planes 
fail  and  the  column  of  disks  rotates  by  rolling  over  the  bottom  of  the  assembly. 
This  is  simulated  by  the  mechanism  framework  shown  in  Fig.  3.11b.  Appropriate 
boundary  conditions  are  simulated  through  rollers,  hinges  and  rigid  boimdary 
trusses. 

The  incremental  shear  strain  applied  on  the  lattice  equivalent  of  the 
particulate  assembly  in  Fig.  3.11a  is  0.00001,  the  reason  for  ihis  limiting  value  was 
discussed  earlier.  However,  since  the  mechanism  approach  is  a  method  to 
simulate  the  predominant  deformation  due  to  sliding  and  rotation  of  particles,  the 
incremental  shear  strain  on  the  mechanism  framework  (Fig.  3.11b)  is  0.001  which 
is  100  times  0.00001.  At  a  rate  of  the  above  incremental  shear  strain,  the  niimber 
of  increments  required  to  reach  a  shear  strain  of  0.20  will  be  0.20/0.001  =  200.  It 
was  shown  in  the  previous  section  that  the  number  of  increments  required  to 
reach  a  shear  strain  of  0.20  is  20,000  using  the  single  particle  approach.  Therefore, 
it  can  be  seen  that  the  mechanism  approach  which  takes  only  200  increments  is 
computationally  very  efficient. 

With  the  centroids  of  the  particles  updated  through  the  mechanism  action, 
the  lattice  shown  in  Fig.  3.1  is  reformed  in  the  new  position.  In  this  new  position, 
the  axial  forces  in  the  bars  prior  to  the  mechanism  are  input  as  initial  bar  forces 
and  an  initial  load  analysis  is  carried  out  for  equilibrium.  The  new  shear  stress 
ratios  at  the  contacts  are  then  computed  and  incremental  load  is  applied.  The 
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criterion  for  mechanism  is  again  checked  and  the  process  is  repeated. 

3.3  Numerical  implementation  of  the  lattice  type  model 

3.3.1  General 

A  computer  program  was  written  to  implement  the  concepts  of  the  lattice 
type  model  in  section  3.2.  Two  different  data  set  are  prepared:  truss  data  and 
particles  data  which  is  referred  to  as  disks  data  here.  Truss  data  consist  of  nodal 
coordinate  data,  connectivity  data,  boimdary  condition  data,  geometric  property 
data,  material  property  data  and  load  data.  Disks  data  consist  of  the  coordinates 
of  the  centroids  of  disks,  number  of  contacts  around  a  disk,  number  of  bars  in  a 
disk,  node  numbers  of  the  contact  nodes,  member  numbers  linked  to  the  contact 
nodes,  back  nodes  of  the  bars  linked  to  the  contact  nodes,  disk  numbers  linked 
to  the  contact  nodes  and  the  node  and  disk  numbers  linked  to  the  boundaries. 
The  authors  has  prepared  computer  programs  to  generate  the  data  for  regular 
arrays  of  circular  disks  arranged  as  very  loose  and  very  dense  packings  (Fig. 
3.12).  The  basic  data  required  to  prepare  the  above  data  are  the  radii  of  the 
particles  and  the  X-Y  coordinates  of  the  centroid  of  the  particles. 


(a)  Loose  packing  (96  disks)  (b)  Dense  packing  (104  disks) 


Fig.  3.12  Loose  and  Dense  packings 

With  the  above  data,  an  assembly  of  particles  transformed  into  a  truss  is 
analyzed  for  a  given  load  and  prescribed  loading  condition  (biaxial,  simple  shear 
etc.)  along  the  procedures  listed  in  the  flow  chart  shown  in  Fig.  3.13. 
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Fig.  3.13  Flow  chart  of  procedures  in  the  Lattice  Type  Model 
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3.3.2  Techniques  in  numerical  simulation  -  A  brief  outline 

Several  numerical  problems  were  dealt  with  during  the  course  of  this 
research  and  these  are  listed  here  as  general  computational  information  to  benefit 
future  researchers  in  lattice  modeling. 

The  boundaries  for  simulating  non-ideal  simple  shearing  of  an  assembly 
of  disks  shown  in  Fig.  3.11a  is  simulated  through  a  top  rigid  truss  supported  by 
two  horizontal  rollers  at  the  two  bottom  ends  of  the  rigid  truss  to  allow  the  free 
vertical  movement  of  the  top  boundary.  The  elastic  modulus  of  each  bar  of  this 
truss  is  set  equal  to  the  elastic  modulus  of  the  particles.  The  cross-sectional  areas 
of  the  bars  of  this  truss  are  set  equal  to  10*®  times  the  maximum  of  all  the  cross- 
sectional  areas  of  the  bars  simulating  the  particulate  assembly.  In  the  lattice  type 
model  (LTM),  the  axial  forces  are  redistributed  for  equilibrium  after  updating  the 
coordinates  of  the  nodes  of  the  truss.  At  this  stage,  the  axial  forces  in  the  rigid 
truss  are  not  redistributed  as  these  will  have  considerable  influence  on  the 
particulate  assembly  due  to  it's  large  stiffness.  By  setting  the  two  ends  of  the 
rigid  truss  on  rollers,  it  was  foimd  that  the  top  truss  does  not  remain  horizontal 
and  deflects  as  shown  in  Fig.  3.14  when  a  shear  strain  is  applied.  Originally,  it 
was  brought  to  horizontal  position  by  applying  a  vertical  displacement  for  the 
two  ends  of  the  truss  (Fig.  3.14).  Thus,  a  truss  analysis  was  to  be  carried  out  to 
bring  the  top  truss  to  a  horizontal  position  after  every  analysis  for  incremental 
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load,  redistribution  of  additional  shear  forces  etc. 


This  was  a  time  consuming  process  and  thus  the  boundary  conditions  were 
changed  to  fix  every  node  of  the  top  rigid  truss  in  the  horizontal  direction.  With 
these  boundary  conditions,  it  was  found  that  the  vertical  displacements  of  the  top 
truss  nodes  were  accurate  to  six  decimal  places  and  thus  the  top  truss  did  not 
deflect. 

Frictionless  boundaries  simulated  by  rollers  as  shown  in  Fig.  3.11a  has 
boundary  effects  and  thus  the  planes  adjacent  to  the  vertical  boundaries  never 
fail.  Therefore,  the  behavior  of  a  particulate  assembly  will  not  be  very  accurate 
in  regions  with  at  least  one  to  two  particle  diameters  away  from  the  vertical 
boundaries.  This  leads  to  particles  not  separating  (contact  forces  becoming  zero) 
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even  after  the  distance  between  the  particles  exceeds  the  sum  of  the  radii  of  the 
particles  in  the  new  position  after  updating  the  coordinates  through  the 
mechanism  framework.  To  rectify  this  problem,  the  contacts  are  split  and  the 
axial  forces  are  redistributed  when  the  distance  between  two  particles  exceeds  the 
sum  of  their  radii  plus  a  tolerance  of  2%  of  the  sum  of  radii. 

When  the  axial  forces  in  the  bars  become  tensile  forces,  a  very  low  cross- 
sectional  area  is  set  for  those  bars  (AWBAR  =  1 0'®  times  the  maximum  of  all  the 
cross-sectional  areas  in  the  particulate  assembly).  When  these  weak  bars  get 
linked  to  the  top  rigid  truss  at  the  boundary,  there  is  a  difference  in  stiffness 
between  these  bars  and  the  bars  of  the  rigid  truss  with  an  order  of  magnitude  of 
10'^  This  would  be  the  maximum  order  of  difference  that  can  be  obtained  with 
computations  accurate  to  a  double  precision.  However  even  with  a  bar  area  of 
AWBAR  for  the  weak  bars,  axial  forces  upto  a  maximum  of  10"^  times  the 
maximum  of  the  axial  forces  in  the  bars  continues  to  develop  in  the  weak  bars. 
Therefore,  these  bar  forces  are  set  to  zero  and  an  initial  stress  analysis  is  carried 
out  for  equilibrium. 

As  particles  lose  contact,  stress-free  zones  are  developed,  where  particles 
remain  in  contact  with  neighboring  particles  or  boundary  with  just  one  or  two 
contacts.  If  a  particle  has  at  least  two  contacts  linked  to  the  global  truss  shown 
in  Fig.  3.1,  it  is  sufficient  to  keep  the  stiffness  matrix  non-singular.  However, 
even  if  only  two  contacts  of  a  particle  are  in  contact  with  neighboring  particles. 
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there  are  situations  of  the  stiffness  matrix  becoming  singular.  One  such  situation 
is  shown  in  Fig.  3.15.  From  this  figure  it  can  be  seen  that  three  particles  can 
group  together  and 


^Group  of  three  particles  separating 
from  the  main  particle  assembly 


Fig.  3.15  Particles  separating  from  the  main 
particulate  assembly 

can  get  separated  from  the  global  truss  or  the  stressed  part  of  the  particulate 
assembly  resulting  in  the  singularity  of  the  [K]  matrbc.  Therefore,  a  condition  is 
stipulated  that  at  least  three  nodes  of  a  particle  should  be  in  contact  witii 
neighboring  particles,  and  the  areas  of  the  bars  linked  to  these  nodes  are  set  equal 
to  AWBAR. 

Because  of  the  above  mentioned  constraints  and  boundary  effects,  the 
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accuraq/^  of  computations  decreases  and  the  errors  get  accumulated  as  a 
particulate  assembly  undergoes  deformation.  Therefore,  contact  forces  do  not 
reduce  to  exactly  zero,  but  only  to  a  very  low  value.  Hence,  a  critical  value  is  set 
for  the  contact  force  in  order  to  split  a  node.  The  critical  value  of  contact  force 
at  which  a  contact  is  split  is  increased  as  the  material  undergoes  deformation. 
The  critical  value  is  increased  as  constraint  ratio  approaches  one.  The  variation 
in  the  critical  value  is  due  to  the  decrease  in  the  accuracy  of  computations  during 
a  deformation.  The  initial  value  is  set  as  0.10%  of  tihe  maximum  contact  force  and 
is  increased  to  2.0%  of  the  maximum  contact  force  when  constraint  ratio  equals 
one. 

It  was  stated  that  a  mechanism  is  identified  when  the  increment  in  load  is 
equal  to  zero.  However,  the  numerical  tests  reveal  that  it  does  not  always 
become  zero  but  gives  a  low  positive  value.  This  is  assumed  to  be  due  to  the 
presense  of  the  existence  of  weak  bars  at  forced  contacts  and  boundary  effects 
explained  earlier.  Hence  a  mechanism  is  identified  if  the  increment  in  load  is  less 
than  a  minimum  value.  This  value  is  equal  0.50%  of  the  applied  load  (axial  load 
in  a  biaxial  test  and  shear  force  in  a  simple  shear  test)  from  zero  strain  to  the 
incremental  strain  in  the  initial  tangent  of  the  stress-strain  curve  for  the  material. 

In  the  mechanism  framework  shown  in  Fig.  3.11b,  it  was  stated  that  weak 
bars  are  introduced  by  linking  particle  centroids  if  the  particles  have  a  sliding 
contact.  A  contact  is  a  stick  and  slip  contact  if  the  distance  between  two  particles 
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is  less  than  the  sum  of  the  radii  of  the  particles  plus  Ax  (0.01  mm  =  0.2%  of  the 
sum  of  radii).  Particles  with  such  stick  and  slip  contacts  are  linked  by  weak  bars 
in  the  mechanism  framework. 

When  an  initial  stress  analysis  is  carried  out  after  updating  the  coordinates 
or  updating  the  cross-sectional  areas  of  the  bars,  the  truss  in  Fig.  3.1  is  analyzed 
for  imbalanced  forces  coming  at  the  nodes.  It  has  been  found  that  even  if  a  very 
small  unbalanced  force  (say  O.OOOIN)  acts  at  a  node  where  there  are  only  two 
strong  bars  (bar  area  >  100.00  x  AWBAR),  it  causes  enormous  displacement  and 
alters  the  results  drastically.  Therefore  a  condition  is  set  that  if  there  are  only  two 
active  bars  (bar  area  >  100.00  x  AWBAR)  linked  to  a  node,  these  bar  areas  should 
be  greater  than  1 0^  times  AWBAR  and  also  that  the  minimum  of  the  angles  made 
by  the  bar  with  X-axis  and  Y-axis  should  be  greater  than  10  degrees.  Otherwise, 
the  unbalanced  loads  at  these  nodes  are  set  to  zero  and  it  has  been  found  that  this 
leads  to  an  accuracy  loss  in  the  range  of  ±  2%. 

3.4  Lattice  modeling  with  centroid  of  particles 

The  authors  began  the  work  by  simulating  the  lattice  by  linking  the  particle 
contacts  with  the  particle  centroids  as  the  nodes  (Fig.  3.16).  Even  though  the 
results  were  in  general  agreement  with  the  behavior  of  granular  media,  it  was 
found  that  the  method  was  not  sufficiently  general  and  can  lead  to  sensuous 


errors.  Therefore,  the  modeling  approach  was  changed  to  a  lattice  formed  by 
linking  only  the  particle  contacts,  the  results  of  which  will  be  the  subject  matter 
of  this  research.  Considerable  work  was  done  using  the  modeling  approach 
shown 


(a)  Particle  arrangement 


(b)  Lattice  type  simulation 


Fig.  3.16  Lattice  type  simulation  with  particle  centroids 


in  Fig.  3.16.  The  results  and  methods  adopted  in  this  approach  and  it's 
shortcomings  are  presented  in  Budhu  et.al  (1994)  and  Ramakrishnan  (1996). 
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3.5  Differences  in  approach  between  Lattice  Type  Model  and 
Discrete  Element  Models 

A  distinct  feature  of  the  lattice  type  model  is  that  the  particle  itself  is 
replaced  by  a  lattice  and  thus  the  (elastic)  behavior  of  the  particle  is  retained 
within  the  microstructure  of  the  particle.  In  the  current  numerical  models  for 
particulate  media  (DEM  and  DMA),  the  elastic  behavior  of  the  particles  are 
replaced  by  virtual  springs  at  contacts,  the  disks  are  translated  and  rotated 
rigidly.  The  lateral  forces  developed  due  to  the  Poisson  effect  of  the  particles  are 
not  neglected  in  the  lattice  type  model. 

In  the  DEM  and  DMA,  the  stiffnesses  of  the  springs  are  a  function  of 
elastic  modulii  and  radii  of  the  disks  in  contact.  If  two  particles  of  different  radii 
are  in  contact,  then  the  contact  stiffnesses  of  the  springs  are  evaluated  using  an 
average  radius.  However,  in  the  lattice  type  model  the  stiffnesses  of  the  bars  of 
the  lattice  are  a  function  of  the  radius  of  the  particle  forming  the  lattice  and  there 
is  no  averaging  procedure.  In  the  available  literature  on  DEM,  there  is  no 
information  as  how  to  compute  the  contact  stiffnesses  if  two  particles  of  different 
elastic  modulii  are  in  contact.  In  the  lattice  type  model  the  elastic  modulus  of  the 
bars  of  the  lattice  being  equal  to  the  elastic  modulus  of  the  particles,  the  stiffness 
of  the  lattice  is  proportioned  to  the  stiffness  of  the  particles.  Also  in  the  lattice 
type  model,  as  the  contacts  are  modeled  as  joints  in  a  lattice,  the  interaction  at 
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contacts  will  automatically  get  related  to  the  relative  stiffriesses  of  the  disks  in 
contact.  Thus,  the  lattice  tjqje  model  has  significant  potential  to  study  the  effects 
of  weak  particles  (particles  with  low  elastic  modulus)  randomly  placed  within  a 
particulate  medium,  on  it's  macro-mechanical  behavior. 

A  notable  feature  of  the  lattice  type  model  is  that  the  relative  sliding, 
rotation  and  rolling  of  particles  are  all  handled  in  one  step  through  a  global 
approach  using  the  mechanism  truss.  Also,  the  development  of  a  mechanism  is 
automatically  highlighted  by  the  global  stiffness  matrix  (due  to  the  reformation 
of  the  bar  areas)  through  the  increment  in  load  becoming  zero. 

Although  the  potential  of  the  lattice  type  model  to  carry  out  study  the 
behavior  of  particulate  media  with  arbitrary  particle  shapes  is  not  yet  carried  out, 
it  can  be  stated  that  the  method  has  capabilities  for  this.  The  lattice  formation  of 
a  random  packing  can  be  carried  out  by  taking  an  image  of  the  particle  assembly. 
The  bar  areas  are  to  be  evolved  by  taking  the  nearest  ellipse  fitting  the  particle 
shape.  Also,  the  model  has  capabilities  to  study  the  effect  of  particle  fracture. 
Application  of  the  lattice  type  model  for  random  packings  with  random  particle 
shapes  and  for  study  of  particle  fracture  shall  be  a  scope  for  future  research. 


CHAPTER  4 


APPLICATION  OF  LATTICE  TYPE  MODEL 
TO  A  TWO  DIMENSIONAL  ARRAY  OF  DISKS 

4.1  General 

In  this  chapter,  the  lattice  type  model  is  applied  to  two  dimensional 
assemblies  of  circular  disks.  The  disks  of  diameter  2.54  mm  are  arranged  into  a 
loosest  packing  of  96  particles  and  a  densest  packing  of  104  particles  (Fig.  3.12). 
The  material  of  the  spheres  is  considered  to  be  quartz.  The  packings  are 
subjected  to  an  initial  vertical  load  followed  by  a  simple  shear  strain  loading  (Fig. 
4.1). 

4.2  Material  properties 

The  properties  of  the  material  quartz  considered  herein  are  (Mantell,  1958) 
Elastic  modulus  (E)  =  50.0  GPa. 

Poisson's  ratio  (u)  =  0.20 

Shear  modulus  (G)  =  20.8  GPa 

Friction  coefficient  (p)  =  0.50 
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vertical  load 


Fig.  4.1  Simulation  of  simple  shear 


4.3  Geometric  properties 


The  geometric  properties  required  for  the  analysis  of  a  truss  are  the  cross- 
sectional  areas  and  lengths  of  the  bars.  The  lengths  of  bars  are  determined  from 
the  X  and  Y  coordinates  of  the  contact  nodes.  The  initial  cross-sectional  areas  of 
the  bars  are  obtained  by  procedures  described  in  Chapter.  3. 
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4.4  Boundary  conditions 

For  both  the  loose  and  dense  packings  the  following  boundary  conditions 
are  assumed.  (1)  The  bottom  bormdary  is  fixed  and  can  mobilize  friction  and 
thus  the  bottom  nodes  are  taken  to  be  on  hinges.  (2)  The  vertical  boimdaries  are 
frictionless  and  are  simulated  by  rollers.  (3)  The  top  boundary  is  simulated  by 
a  very  rigid  truss  to  represent  a  load  platen.  The  nodes  of  the  rigid  truss  are 
restrained  in  the  X-direction  and  are  free  in  the  Y-direction.  Thus,  the  top  truss 
is  allowed  to  move  up  or  down,  simulating  the  so-called  non-ideal  simple  shear 
conditions  (Budhu,  1988),  i.e.,  allowing  a  volume  change  of  the  material. 

4.5  Loading  condition  and  Lattice/ truss  analysis 

The  loading  is  comprised  of  a  vertical  load  followed  by  the  superposition 
of  a  simple  shear  displacement  (forward  shear)  as  shown  in  Fig.  4.1.  The  total 
applied  vertical  load  is  268.00  N.  The  vertical  load  is  applied  as  nodal  loads  of 
22.16  N  at  the  top  level  nodes  of  the  rigid  truss  (Fig.  4.1).  The  simple  shear 
straining  is  simulated  by  applying  a  triangular  variation  of  displacements  for  the 
nodes  attached  to  the  vertical  boundaries  as  shown  in  Fig.  4.1. 

The  particulate  packings,  simulated  as  a  lattice/truss  as  shown  in  Fig.  3.1b, 
are  analyzed  for  the  loading  condition  along  the  steps  given  in  the  flow  chart 
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shown  in  Chapter.  3  (Fig.  3.13).  The  results  of  the  analysis  of  the  loose  and  dense 
packings  are  presented  separately  in  Sections  4.7  and  4.8.  The  results  are 
presented  in  relation  to  the  state  of  stress  in  simple  shear  and  an  average  stress 
tensor  which  are  briefly  presented  in  the  next  section. 

4.6  State  of  stress  in  simple  shear  and  average  stresses 

The  ideal  system  of  stresses  on  an  element,  rectangular  in  cross-section 
subjected  to  simple  shear  strain  is  shown  in  Fig.  4.2.  In  this  figure,  the  following 
notations  are  used  for  the  stresses.  0^^^^  -  Normal  stress  on  vertical  planes,  - 

Normal  stress  on  horizontal  planes,  -  Shear  stress  on  vertical  planes,  - 

Shear  stress  on  horizontal  planes,  O22  -  Lateral  stress.  In  general,  it  is  difficult  to 
reproduce  these  stresses  in  the  laboratory  tests,  in  particular  the  complementary 
shear  stresses  and  Consequently,  the  distributions  of  stresses  and  strains 
near  the  vertical  boundaries  of  laboratory  samples  are  non-uniform.  However, 
the  state  of  stress  in  the  central  zone  of  particulate  packings  (Examples:  Figs.  4.3a 
and  4.3b)  are  expected  to  be  uniform.  The  state  of  stress  in  an  infinitesimally 
small  element  at  the  center  of  the  central  core  can  be  assumed  to  represent  the 
state  of  stress  in  the  central  core.  The  components  of  this  stress  tensor  in  the 
central  core  can  be  evaluated  from  the  magnitude  and  location  of  the  discrete 
contact  forces  acting  at  the  boimdaries  of  the  central  core. 
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The  concept  of  evaluating  an  average  stress  tensor  to  determine  the  average 
constitutive  behavior  of  a  certain  region  of  a  material  was  addressed  by  Hill 
(1963).  The  procedure  is  to  convert  a  discrete  set  of  contact  forces  to  an 


Fig.  4.2  Stress  for  Ideal  simple  shear 


Fig.  4.3  Central  core  in  (a)  Loose  packing  (b)  Dense  packing 
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equivalent  state  of  stress  in  a  continuum.  Based  on  Hill's  averaging  principle,  the 
average  stresses  in  a  two-dimensional  element  O^,  T^„  and  at  the  center 
of  tile  particulate  packing  shown  in  Fig.  4.1  can  be  computed.  Hill's  concept  was 
also  applied  by  Drescher  and  De  Josselin  de  Jong  (1972)  for  the  evaluation  of 
average  stresses  in  assemblies  of  disks  using  the  contact  forces  determined  from 
photoelastic  experiments.  The  average  stress  tensor  O^.  at  the  center  of  a  volume 
of  material  shown  in  Fig.  4.4  with  a  Cartesian  coordinate  system  1,  2  is  defined 
as: 


=  p  /,  "s  ^ 

where  0,y  is  the  actual  stress  state  at  each  point  of  the  body  and  V  is  its  volume. 
The  average  stress  tensor  is  defined  for  an  infinitesimally  small  element  at  the 
center  of  the  control  volume  and  it  can  be  transformed  into  an  expression 
containing  the  forces  acting  at  the  boundary  of  the  control  volume  as  follows. 

Oj-  =  ^ik  °ik  =  hk  ^kj  (4.2) 

where  5;^  is  the  Kronecker  delta,  x,.  are  the  coordinate  directions  (i  =  1,  2,  3  for 
X,  y,  z  directions)  and  the  comma  indicates  partial  differentiation  with  respect  to 
the  coordinates.  Substituting  Equation.  (4.2)  into  Equation  (4.1),  we  get 


Fig.  4.4  General  distribution  of  forces  and  Tj  around  surface  S  of  volume  V 
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a,j  dV 


(4.3) 


The  stress  state  satisfies  the  equilibrium  equations 


■'*/.  j 


=  0 


(4.4) 


Using  this  condition.  Equation.  (4.3)  can  be  expressed  as 


(4.5) 


Using  Gauss's  divergence  theorem.  Equation.  (4.5)  can  be  written  as 
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(4.6) 


where  Hj  are  the  vector  of  direction  cosines  of  the  outward  normal  at  coordinates  Xj 
on  the  surface  and  Pj  are  the  vector  of  stress  components.  Thus,  it  can  be  seen 
from  Equation.  (4.6)  that  the  average  stress  tensor  can  be  obtained  by  integrating 
the  products  Pj  Xj  over  the  boundary  of  the  volume  V.  For  a  set  of  discrete  forces 
Tj  over  the  surface  (Fig.  4.4),  Equation.  (4.6)  becomes 


a  =  —  Y  T 


m 


(4.7) 


The  average  stresses  at  the  center  of  the  assembly  of  particles  shown  in  Fig. 
4.1  are  foimd  using  Equation.  (4.7)  from  the  normal  and  shear  forces  at  the 
contacts  between  the  particles  and  the  boundaries.  In  the  absence  of  body 
couples  in  a  two  dimensional  element,  =  T^.  The  shear  stress  ratios  on  the 
vertical  and  horizontal  planes  are  defined  as  and  respectively. 

From  the  components  of  tiie  average  stress  tensor,  the  principal  stresses 
and  O33  and  the  angle  i|;  made  by  the  major  principal  stress  plane  with  the 
positive  direction  of  X-axis  (Fig.  4.5)  are  determined  from 
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(4.8) 
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4.7  Loose  packing,  results  and  discussions 

4.7.1  Shear  stress  ratios  and  contact  forces 

The  shear  stress  ratios  at  the  contact  nodes  for  the  loose  packing  are  shown 
in  Figs.  4.6a  and  4.6b  at  small  shear  strains  (y)  of  0.001%,  0.002%  and  in  Figs  4.6c 
and  4.6d  at  large  shear  strains  of  30%  and  40%.  The  above  shear  strains  are 
chosen  in  order  to  show  the  large  variation  in  the  magnitudes  of  the  shear  stress 
ratios  during  the  shearing  process.  Relative  sliding  of  particles  occurs  at  a  contact 
when  the  shear  stress  ratio  at  that  contact  exceeds  the  coefficient  of  friction. 
When  there  is  a  continuous  line/curve  of  sliding  contacts,  a  slippage /sliding 
plane  develops.  A  comparison  of  Fig.  4.6a  and  Fig.  4.6b  shows  that  slippage  is 
initiated  on  vertical  planes  and  the  deformation  pattern  quickly  resembles  the 
tilting  and  frictional  sliding  of  a  column  of  particles.  This  deformation  pattern  or 
mechanism  is  then  simulated  by  the  "mechanism  framework"  explained  in 
Chapter.  3.  For  examples,  the  mechanism  frameworks  at  shear  strains  of  0.002% 
and  30%  are  shown  in  Figs.  4.7a  and  4.7b.  The  deformation  is  one  of  a  framed 
pin  jointed  structure  subjected  to  lateral  loads. 

The  change  in  the  magnitude  of  normal  forces  at  contacts  are  shown  in 
Figs.  4.8a,  4.8b  and  4.8c.  The  change  in  normal  forces  on  the  horizontal  contact 
planes  is  only  marginal  while  there  is  a  substantial  increase  in  the  normal  forces 
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Fig.  4.6a  Shear  stress  ratios  at  contacts  for  loose  packing  at  y  =  0.001  % 
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4.6b  Shear  stress  ratios  at  contacts  for  loose  packing  at  y  =  0.002  % 


Fig.  4.6d  Shear  stress  ratios  at  contacts  for  loose  packing  at  y  =  40  % 
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(a) 


(b) 


Fig.  4.7a, b  Mechanism  frameworks  for  loose  packing  at 
(a)  Y  =  0.002%  (b)  y  =  30% 
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Fig.  4.8b  Normal  forces  at  contacts  for  loose  packing  at  y  =  10  % 
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Fig.  4.8c  Normal  forces  at  contacts  for  loose  packing  at  y  =  40  % 
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on  the  vertical  contact  planes.  Since  the  contacts  on  vertical  planes  continue  to 
slide  and  the  normal  forces  increase,  the  shear  forces  at  these  contacts  increase 
considerably.  For  equilibrium  of  the  particles,  the  shear  forces  at  the  contacts  on 
horizontal  planes  should  also  increase.  This  can  be  seen  from  the  substantial 
increase  in  shear  stress  ratios  at  the  contacts  on  horizontal  planes  (Fig.  4.6c)  at  a 
large  shear  strain  of  30%,  while  there  is  only  a  marginal  decrease  in  the  normal 
forces.  As  the  material  tends  to  move  towards  the  closest  packing  configuration, 
the  shear  stress  ratios  at  contacts  on  vertical  planes  continue  to  be  sliding  ones 
until  a  shear  strain  of  about  30%  while  there  is  a  considerable  increase  in  the 
shear  stress  ratios  on  horizontal  planes  due  to  the  increased  shear  forces  at 
contacts  on  these  planes.  At  a  shear  strain  of  40%  (Fig.  4.6d),  the  shear  stress 
ratios  at  contacts  on  vertical  planes  decrease  and  they  now  become  non-sliding 
contacts.  At  this  shear  strain,  it  can  be  seen  that  continuity  of  sliding  contacts 
prevails  neither  along  horizontal  planes  nor  along  vertical  planes. 

4.7.2  Contact  forces  and  fabric  changes 

The  variation  of  the  magnitude  and  orientation  of  contact  forces  are  shown 
using  vector  plot  of  contact  forces  within  the  particulate  assemblies  (Fig.  4.9a). 
The  thickness  of  lines  in  these  diagrams  indicates  the  magnitude  of  the  contact 
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Fig.  4.9a,b  Vector  plot  of  contact  forces  for  loose  packing  under 
(a)  Vertical  load  (b)  y  =  0.002  % 
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Fig.  4.9c,d  Vector  plot  of  contact  forces  for  loose  packing  at 
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Fig.  4.9e  Vector  plot  of  contact  forces  for  loose  packing  at  y  =  50% 
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forces  relative  to  the  maximum  cor\tact  force  in  the  particulate  assembly.  These 
diagrams  show  the  flow  of  contact  forces  and  load  paths  as  well  as  the  changes 
in  fabric  or  packing  structure. 

The  vector  plot  of  contact  forces  are  shown  in  Fig.  4.9a  imder  vertical  load 
and  in  Figs.  4.9b,  4.9c,  4.9d  and  4.9e  at  shear  strains  of  0.002%,  10%,  30%  and  50% 
respectively.  Obliquity  of  contact  force  is  defined  as  the  angle  made  by  the 
orientation  of  the  contact  force  with  the  contact  normal  at  the  contact.  The 
contact  normal  for  two  particles  of  circular  cross-section  is  the  direction  linking 
the  centroid  of  the  particles.  The  obliquity  of  contact  forces  on  vertical  planes 
imder  vertical  load  which  were  initially  horizontal  (Fig.  4.9a)  is  inclined  at  an 
angle  of  26.5°  (tan”^0.5)  at  a  shear  strain  of  0.002%  (Fig.  4.9b).  At  a  shear  strain 
of  10%  (Fig.  4.9c),  the  magnitude  and  orientation  of  contact  forces  in  the  middle 
ten  coluiims  of  particles  are  imiform  and  non-imiformity  can  be  seen  only  along 
the  column  of  particles  adjacent  to  the  vertical  boundaries.  However,  at  a  larger 
shear  strain  of  30%  (Fig.  4.9d),  contact  forces  of  larger  magnitude  can  be  seen  in 
the  triangular  zones  to  the  right  of  line  D2-D2  and  to  the  left  of  line  Dl-Dl.  This 
non-uniformity  is  likely  to  be  caused  by  the  interaction  between  the  right 
boundary  and  top  boimdary  and  between  the  left  boimdary  and  bottom  boundary 
through  the  triangular  zones  shown  in  Fig.  4.9d. 

The  material  tended  to  reach  the  closest  packing  configuration  at  a  shear 
strain  of  50%  (Fig.  4.9e).  At  this  shear  strain,  it  can  be  seen  that  new  contacts 
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have  been  established  along  columns  of  particles  parallel  to  line  Dl-Dl.  The 
shearing  could  not  be  carried  out  further,  as  the  scope  of  the  computer  code  in 
this  research  does  not  incorporate  the  formation  of  new  contact  nodes  and 
elements  to  develop  a  revised  lattice. 

From  the  vector  plot  of  contact  forces,  it  can  be  seen  that  over  60%  of  the 
assembly  in  the  middle  is  under  tmiform  loads.  At  the  ends,  the  loads  are 
significantly  higher  or  lower  (depending  on  the  direction  of  shear)  than  the  loads 
in  the  middle  60%  of  the  assembly.  This  is  expected  since  the  vertical  sides  of  the 
particulate  assembly  are  assumed  to  be  frictionless  so  that  no  complementary 
shear  stress  is  developed  there.  The  contact  force  distributions  from  the  LTM 
results  are  similar  to  boimdary  stress  distributior\s  in  non-ideal  simple  shear  using 
elastic  analysis  of  a  continuum  discussed  by  Roscoe  (1953). 

At  higher  shear  strains,  the  contact  force  flow  pattern  resembles  a 
staircase /steplike  pattern.  The  contact  forces  on  vertical  planes  increase  and 
become  almost  equal  to  the  contact  forces  on  the  horizontal  planes  as  the  material 
tends  to  reach  the  closest  packing  at  50%.  This  indicates  that  as  the  material 
tends  to  become  a  dense  packing,  the  stress  conditions  reach  a  hydrostatic  state. 
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4.8  Dense  packing,  results  and  discussions 

4.8.1  General 

For  ease  of  following  the  deformation  mechanism  in  the  dense  assembly 
of  disks,  certain  symbols  are  used  for  the  orientation  of  rows  and  columns  of 
particles  and  these  are  shown  in  Fig.  4.10.  Two  lines  AA  and  BB  along  the  center 
line  of  diagonal  particles  are  selected  with  line  AA  representing  the  short 
diagonal  (compression)  and  line  BB  representing  the  long  diagonal  (tension) 
during  simple  shear.  The  horizontal  row  of  particles  are  labelled  as  shown  in  Fig. 
4.10. 

4.8.2  Shear  stress  ratios,  contact  forces  and  fabric  changes 

The  shear  stress  ratios  at  the  contact  nodes  are  shown  in  Fig.  4.11a  under 
vertical  load  and  in  Figs.  4.11b,  4.11c,  4.11d  and  4.11e  at  shear  strains  of  0.01%, 
0.03%,  0.04%  and  10%  respectively.  The  normal  component  of  the  contact  forces 
at  contact  nodes  are  shown  in  Fig.  4.12a  under  the  vertical  load  and  in  Figs.  4.12b, 
4.12c,  4.12d  and  4.12e  at  shear  strains  of  0.01%,  0.03%,  0.04%  and  10%.  The 
dense  packing  arrangement  in  Fig.  4.1  is  the  densest  packing  that  is  physically 
possible  in  an  assembly  of  disks  of  imiform  size.  Therefore,  the  material 
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Fig.  4.11b  Shear  stress  ratios  at  contacts  for  dense  packing  at  y  =  0.01  % 
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Fig.  4.11d  Shear  stress  ratios  at  contacts  for  dense  packing  at  y  =  0-04  % 
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Fig.  4.12a  Normal  forces  at  contacts  for  dense  packing  imder  vertical  load 
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Fig.  4.12b  Normal  forces  at  contacts  for  dense  packing  at  y  =  0.01% 
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Fig.  4.12d  Normal  forces  at  contacts  for  dense  packing  at  y  =  0.04  % 
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Fig.  4.12e  Normal  forces  at  contacts  for  dense  packing  at  y  =  10  % 
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subjected  to  a  simple  shear  strain  has  to  dilate  in  order  to  imdergo  a 
deformation.  The  shear  stress  ratios  under  vertical  load  in  Fig.  4.11a  shows  that 
(n-2)  contacts  are  sliding  in  most  of  the  particles  where  'n'  is  the  number  of 
contacts  around  a  particle.  However,  there  is  no  mechanism  at  this  stage  as  the 
material  can  carry  deviatoric  stress  as  will  be  seen  in  the  stress-strain  behavior 
later.  Also,  the  constraint  ratio  computed  under  such  condition  is  1.7  which  is  far 
greater  than  one,  a  condition  for  the  presense  of  a  mechanism.  An  examination 
of  Figs.  4.11a,  4.11b  and  4.11c  shows  the  progress  of  shear  stress  ratios  leading  to 
sliding  contacts  at  the  contacts  on  vertical  planes.  A  comparison  of  Figs.  4.12a 
and  4.12b  shows  that  the  first  step  is  the  loss  of  contacts  aligned  along  the  lines 
Dl-Dl,  D2-D2,  D3-D3  and  D4-D4.  It  can  also  be  seen  that  the  normal  forces  along 
the  line  of  contacts  parallel  to  line  AA  increase  from  11.0  N  to  an  average  value 
of  18.0  N,  while  the  normal  forces  along  the  line  of  contacts  parallel  to  line  BB 
decreases  from  11.0  N  to  an  average  value  of  8.0  N.  In  summary.  Figs.  4.12a, 
4.12b  and  4.12c  shows  the  increase  or  decrease  of  normal  forces  leading  to 
stronger  contacts  or  bondless  nodes. 

Through  the  above  mentioned  transformations,  the  assembly  develops  into 
a  mechanism  at  a  shear  strain  of  0.04%.  Further  deformations  are  simulated 
through  the  mechanism  frameworks  shown  in  Figs.  4.13a,  4.13b,  4.13c  and  4.13d 


at  various  shear  strains. 


•  PflRTiaX  CENTROID 


Fig.  4.13a,b  Mechanism  frameworks  for  dense  packing  at 
(a)  Y  =  0.04  %  (b)  Y  =  1  % 
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Fig.  4.13c,d  Mechanism  frameworks  for  dense  packing  at 
(c)  Y  =  10  %  (d)  Y  =  20  % 
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Fig.  4.14a,b  Vector  plot  of  contact  forces  for  dense  packing  under 
(a)  vertical  load  (b)  y  =  0.01  % 
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Fig.  4.14c,d  Vector  plot  of  contact  forces  for  dense  packing  under 
(c)  Y  =  0.03%  (d)  Y  =  0.04  % 
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Fig.  4.14e,f  Vector  plot  of  contact  forces  for  dense  packing  under 

(e)  Y  =  5  %  (f)  Y  =  10  % 
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Fig.  4.14g,h  Vector  plot  of  contact  forces  for  dense  packing  under 

(g)  Y  =  20  %  (h)  Y  =  30  % 
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The  load  path  or  contact  force  flow  over  a  range  of  shear  strains  is  shown 
in  Fig.  4.14a  through  Fig.  4.14h.  Fig.  4.14a  shows  that  larger  loads  flow 
predominantly  in  the  vertical  direction  under  the  initial  vertical  load.  Under  the 
application  of  a  shear  strain  of  0.01%,  the  path  of  the  larger  loads  immediately 
tilts  towards  line  AA  (Fig.  4.14b)  and  at  a  shear  strain  of  0.04%  the  load  path 
within  the  middle  region  of  the  central  core  is  confined  to  paths  along  line  AA 
(Fig.  4.14d).  The  magnitude  and  orientation  of  the  contact  forces  at  a  large  shear 
strain  of  30%  in  Fig.  4.14h  shows  that  certain  particles  are  not  m  equilibrium. 
This  is  due  to  the  accumulation  of  errors  due  to  numerical  approximations 
discussed  in  Chapter.  3  and  also  due  to  the  mechanism  framework  not  giving  the 
correct  position  of  the  particles  at  very  large  shear  strains. 

An  analysis  of  the  variation  of  shear  stress  ratios  between  the  initial  vertical 
load  (Fig.  4.11a)  and  a  shear  strain  of  0.01%  (Fig.  4.11b)  shows  that  the  line  of 
contacts  parallel  to  line  AA  changes  from  sliding  to  non-sliding  contacts,  which 
could  be  attributed  to  the  increase  in  normal  force  at  these  contacts.  However, 
the  line  of  contacts  parallel  to  line  BB  continues  to  be  sliding  contacts  and  this  is 
caused  by  the  decrease  in  normal  forces  at  these  contacts.  The  line  of  contacts 
along  Cl-Cl,  C2-C2  etc.,  which  were  originally  non-sliding  contacts,  changes  to 
sliding  contacts.  This  is  due  to  the  mobilization  of  shear  forces  along  these 
contacts  as  it  can  be  seen  from  Figs.  4.12a  and  4.12b  that  there  is  no  significant 
change  in  the  normal  forces.  The  tranformations  in  contact  forces  and  the  fabric 
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changes  from  the  initial  vertical  load  to  a  shear  strain  of  0.04%  can  be 
summarized  as  follows: 

•  Stress-free  zones  develop  at  top  right  and  bottom  left  comers  of  the 
particle  assembly  (Fig.  4.12d) 

•  Normal  forces  along  the  line  of  contacts  parallel  to  line  AA  have  doubled 
from  the  initial  values  (Fig.  4.12d). 

•  Line  of  contacts  parallel  to  line  BB  become  bondless  nodes. 

•  Contacts  along  Cl-Cl,  C5-C5  and  at  the  interface  between  the  rows  of 
particles  along  C4-C4  and  D4-D4  and  along  C2-C2  and  Dl-Dl  remain  as 
sliding  contacts. 

•  Considerable  reduction  in  the  number  of  bonded  contacts  when  shear 
strain  reaches  0.04%  (Fig.  4.1 2d). 

A  comparison  of  Fig.  4.11a  with  Fig.  4.11b  shows  the  development  of  void 
spaces  between  particles  along  the  column  of  particles  parallel  to  line  BB. 
However,  the  contacts  along  lines  C2-C2,  C3-C3,  C4-C4  and  Dl-Dl,  D2-D2,  D3- 
D3,  D4-D4  remain  as  "stick  and  slip"  contacts.  The  development  of  sliding  and 
bondless  contacts  shows  that  the  failure  mechanism  is  brought  about  by  the 
rotation  of  particles  in  a  column-like  arrangement  along  with  sliding  interfaces 
with  their  neighboring  columns. 
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4.8.3  Distribution  of  normal  and  shear  forces  at  boundaries 

The  variation  of  normal  forces  and  shear  forces  at  the  bottom  and  top 
boundaries  are  shown  in  Figs.  4.15a  and  4.15b  respectively.  The  distribution  of 
normal  forces  along  the  left  and  right  boxmdaries  is  shown  in  Fig.  4.16  at  various 
shear  strains.  The  following  observations  are  made  from  these  graphs: 

•  The  normal  forces  at  the  bottom  boundary,  decreases  near  the  left 
boundary  and  becomes  zero  at  high  shear  strains,  then  gradually  increases 
to  the  middle  region,  then  remains  constant  in  the  middle  third  region  and 
then  increases  in  the  remaining  one  third  region  towards  the  right 
boundary.  The  variation  at  the  top  boundary  follows  a  similar  pattern 
from  right  to  left. 

•  The  normal  forces  on  the  left  boundary  follow  a  variation  close  to  a  linear 
path  with  zero  at  the  bottom  and  increasing  to  a  maximum  value  at  the 
top.  The  variation  on  the  right  boimdary  follows  a  similar  pattern  from 
top  to  bottom. 

•  The  shear  forces  at  the  bottom  boundary  is  zero  at  the  left  end,  increases 
gradually  to  the  middle  region  and  tend  to  remain  constant. 


X  /  L  (Top  boundary) 
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Fig.  4.15a  Distribution  of  normal  forces  on  the  top  and  bottom  boundaries 

for  the  dense  packing 
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Y  /  H  (Right  boundary) 
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Fig.  4.16  Distribution  of  normal  forces  on  the  vertical  boimdaries 

for  the  dense  packing 
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Fig.  4.18  Distribution  of  stresses  in  simple  shear  fi 

(Roscoe,  1953) 
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The  above  observations  are  qualitatively  in  agreement  with  the  distribution 
of  normal  stresses  existing  at  boundaries  in  simple  shear  tests  on  sands  (Budhu, 
1979)  shown  in  Fig.  4.17. 

The  normal  and  shear  force  distributions  at  the  bottom  boundary  in  Figs. 
4.15a  and  4.15b  indicate  that  even  though  non-uniform  states  of  stress  prevail  at 
the  boxmdaries,  conditions  in  the  middle  of  the  particulate  assembly  are  uniform. 
Roscoe  (1953)  analyzed  an  elastic  material  deforming  in  the  Cambridge  simple 
shear  apparatus  using  Airy's  stress  ftmction.  He  showed  that  even  though  the 
stresses  and  strains  on  the  boundaries  of  the  sample  as  a  whole  may  not  be 
uniform,  conditions  in  the  middle  third  of  the  sample  perpendicular  to  the 
direction  of  shear  can  be  expected  to  be  uniform.  Thus  the  observations  from  the 
lattice  type  model  are  in  agreement  with  results  of  Roscoe's  analysis  shown  in  Fig. 
4.18. 


4.8.4  Fabric  changes,  displacement  fields  and  volume  change  behavior 

From  the  deformation  patterns  shown  in  Figs.  4.14a  through  4.14h,  some 
relationships  can  be  evolved  between  the  fabric  changes,  displacement  fields  and 
volume  change  behavior.  The  displacement  fields  are  shown  in  Figs.  4.19a,  4.19b 
and  4.19c  at  shear  strains  of  0.01%,  5%  and  20%  respectively.  The  variation  of 
volumetric  strain  defined  as  the  change  in  volume  over  the  initial  volume  is 
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Fig.  4.19a,b,c  Displacement  fields  for  dense  packing  at 
(a)  Y  =  0.01%  (b)  Y  =  5%  (c)  y  =  20% 
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Fig.  4.20  Variation  of  shear  stress  ratio  at  bottom  and  volumetric  strain 

for  dense  packing 
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Shear  stress  ratio 


Fig.  4.21  Variation  of  shear  stress  ratios  on  vertical  planes  and 
horizontal  planes  and  void  ratio  for  dense  packing 


Void  ratio  ratio 
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shown  in  Fig.  4.20a  over  a  shear  strain  less  than  0.5%  and  in  Fig.  4.20b  at  large 
shear  strains.  The  variation  of  void  ratio  is  shown  in  Figs.  4.21a  and  4.21b. 

The  volumetric  strain  and  the  void  ratio  curves  show  that  the  material 
dilates  as  expected  for  dense  granular  materials  and  the  overall  behavior  from  the 
two  curves  can  be  summarized  as  follows: 

•  Y  -  10%  -  Linear. 

•  10%  >  Y  ^  30%  -  Convex  curves  with  a  decreasing  slope. 

•  30%  >  Y  ^  50%  -  Volumetric  strain  and  void  ratio  tends  to  remain  constant. 

•  The  void  ratio  of  the  dense  packing  changes  from  an  initial  value  of  0.165 
to  a  final  value  of  0.275  at  shear  strain  50%. 

•  The  final  void  ratio  of  0.275  is  almost  equal  to  0.276,  the  initial  void  ratio 
of  the  loose  packing  which  implies  that  the  dense  packing  has  transformed 
into  a  loose  packing. 

The  major  observations  from  the  displacement  fields  (Figs.  4.19a,b,c)  are 
that  the  displacements  are  predominantly: 

•  in  the  X-direction  at  a  very  low  shear  strain  of  0.01%. 

•  in  the  Y-direction  at  an  intermediate  shear  strain  of  5%. 

•  in  X-direction  at  a  very  large  shear  strain  of  20%. 

The  fact  that  there  is  a  steep  rise  in  the  volumetric  strain  imtil  a  shear 
strain  of  10%  is  explained  by  the  above  variation  of  the  displacement  fields  from 
a  shear  strain  of  0.01%  to  5%.  The  reason  for  the  decrease  in  slope  of  the 
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volumetric  strain  curve  is  explained  by  the  tranformation  in  the  displacement 
fields  from  a  shear  strain  of  5%  to  20%.  Thus,  at  large  shear  strains,  the 
displacement  fields  have  a  larger  component  in  the  X-direction  and  at  very  large 
shear  strains,  the  Y-component  of  the  displacement  fields  become  zero  and  thus 
the  material  ceases  to  dilate. 

The  displacement  fields  which  influence  the  volumetric  strain  are  in  turn 
influenced  by  the  fabric  changes  shown  in  Figs.  4.14a  through  4.14h.  The  fabric 
changes  are  caused  by  the  mechanism  frameworks  shown  in  Figs.  4.13a  through 
4.13d.  A  comparison  of  Figs.  4.13a,b,c  shows  that  the  chain  of  particles  with 
strong  bonds  represented  by  the  strong  bars  rotate  about  the  base  of  the  packing 
resulting  in  a  large  component  of  displacement  in  the  Y-direction  when  the 
mechanism  framework  is  subjected  to  a  deformation.  Flowever,  it  is  to  be  noticed 
that  at  a  large  shear  strain  of  20%  (Fig.  4.13d)  the  above  process  has  not  continued 
and  the  strong  chains  can  be  seen  to  collapse.  This  collapse  mechanism  is 
explained  with  respect  to  Fig.  4.13d. 

The  collapse  is  caused  by  the  preference  of  certain  group  of  particles,  to 
stay  in  contact  or  in  "stick  and  slip  contact",  in  certain  regions.  A  general 
statement  was  made  earlier  that  the  fabric  changes  take  place  in  such  a  way  that 
void  spaces  are  developed  between  particles  along  line  BB.  Flowever,  this  is  not 
true  for  certain  region.  Consider  the  row  of  particles  Dl-Dl  and  C2-C2  in  the 
right  half  of  the  particle  assembly  and  the  row  of  particles  D4-D4  and  C4-C4  in 
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the  left  half  of  the  particle  assembly  (Fig.  4.14g).  Here  a  void  space  is  not 
developed  between  the  particles  aligned  along  line  BB.  Now,  comparing  the 
mechanism  frameworks  in  Figs.  4.13c  and  4.13d,  it  is  clear  that  the  above  row  of 
particles  have  only  translated  in  the  X-direction.  However,  it  is  to  be  noticed  that 
particles  along  the  row  D2-D2  have  rotated  with  respect  to  the  particles  along  row 
C2-C2  (Fig.  4.13d).  Similar  rotations  can  be  seen  in  the  top  left  regions  also. 
Comparing  the  initial  position  of  the  particles  and  their  positions  at  shear  strain 
30%,  the  rotations  of  the  particles  are  as  high  as  30°.  Therefore,  the  conclusion 
of  the  above  discussion  is  that  it  is  this  collapse  mechanism  that  results  in  the 
retardation  of  the  volumetric  strain  at  high  shear  strains. 

4.8.5  Shear  stress  ratio-shear  strain  behavior 

The  shear  stress  ratio-shear  strain  relationship  for  the  dense  packing  at  a 
shear  strain  less  than  0.5%  is  shown  in  Fig.4.20a.  Two  curves  are  shown  for  the 
shear  stress  ratio-shear  strain  relationship.  One  is  the  ratio  of  the  arithmetic 
average  of  the  shear  forces  and  normal  forces  over  the  entire  bottom  boundary 
while  the  other  is  for  the  middle  one  third  of  the  bottom  boundary.  The  reason 
for  this  choice  is  that  the  shear  and  normal  forces  are  not  distributed  uniformly 
at  the  top  and  bottom  boundaries  and  a  uniform  distribution  prevails  only  in  the 
middle  third  of  the  boundary  (Figs.  4.15b  and  4.15b).  Thus,  the  average  forces 
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normally  measured  by  load  transducers  in  laboratory  tests  at  the  bottom 
boimdary  are  likely  to  be  different  from  the  average  forces  at  the  middle  of  the 
bottom  boundary.  This  is  confirmed  in  Figs.  4.20a  where  the  middle  gives  a 
higher  shear  stress  ratio  (about  10%-15%)  than  the  entire  bottom  boimdary.  In 
Fig.  4.20a,  a  "kink"  can  be  seen  in  the  curves  at  a  shear  strain  of  0.02%  which  is 
caused  by  a  large  number  of  contacts  becoming  bondless. 

The  shear  stress  ratio-shear  strain  relationship  at  large  strains  is  shown  in 
Fig.  4.20b.  Using  continuum  analogy,  the  general  response  can  be  summarized 
as  follows: 

•  Y  -  0.02%  -  Linear  elastic. 

•  0.02%  >  Y  ^  0.04%  -  A  quasi-linear  curve. 

•  0.04%  >  y  <  0.50%  -  Perfectly  plastic. 

•  y  >  0.50%  -  Strain  softening. 

The  volumetric  strain  curve  in  Fig.  4.20b  shows  that  in  the  range  of  shear 
strain  from  30%  to  50%,  the  volumetric  strain  tends  to  remain  constant.  This 
shows  that  the  material  shows  a  tendency  to  reach  a  sort  of  "critical  state"  at  30% 
shear  strain.  If  the  material  has  reached  "critical  state"  at  the  shear  strain  of  30%, 
the  stress-strain  curve  is  expected  to  become  horizontal.  Even  though  this 
behavior  is  not  seen,  the  stress-strain  curves  show  a  tendency  to  flatten. 

Hill's  averaging  principle  was  used  to  transform  the  discrete  set  of  contact 
forces  on  the  boundaries  to  an  equivalent  state  of  stress  at  the  center  of  the 
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packing  (Fig.  4.3b).  For  the  purpose  of  comparative  study,  the  average  stresses 
are  computed  from  the  contact  forces  at  the  boundary  of  the,  central  core  and  the 
whole  volume  of  the  packing,  and  both  the  results  are  discussed.  The  shear 
stress  ratios  were  found  for  the  vertical  planes  and  for  the  horizontal 

planes  (T^  /  o^).  Two  levels  of  shear  strains  are  selected  to  illustrate  the 
response  at  small  shear  strains,  y  <  0.5%  (Fig.  4.21a)  and  large  shear  strains,  y  > 
0.5%  (Fig.  4.21b).  The  interpretation  of  the  results  can  be  summarized  as  follows. 

(1)  The  shear  stress  ratio  on  the  vertical  planes  increases  much  more  rapidly 
and  to  significantly  higher  levels  than  on  the  horizontal  planes.  For  the  central 
core,  the  ratio  of  the  shear  stress  ratio  on  the  vertical  planes  to  the  horizontal 
planes  is  about  "3"  at  the  shear  strain  y  =  0.04%,  when  a  mechanism  is  reached 
for  the  first  time. 

(2)  The  central  core  of  the  assembly  gives  shear  stress  ratio  on  the  vertical 
planes  about  10%  greater  than  the  whole  assembly  at  a  shear  strain  of  0.04%  and 
the  difference  increases  at  higher  shear  strains.  This  is  reasonable  because  it  was 
shown  that  the  middle  third  of  the  particle  assembly  is  uniformly  stressed  (Fig. 
4.15b  and  4.15b).  The  difference  for  the  horizontal  planes  is  also  about  10%,  but 
this  difference  continues  imtil  50%  shear  strain.  The  above  behaviors  are  perhaps 
due  to  the  freedom  of  movement  for  particles  in  the  vertical  direction  as 
compared  with  zero  strains  in  the  horizontal  direction. 

(3)  The  shear  stress  ratio  on  vertical  planes  increases  and  reaches  a  peak  value 
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of  1.35  at  a  shear  strain  of  0.04%,  which  is  higher  than  the  coefficient  of  friction 
of  the  material.  Beyond  this  shear  strain,  the  ratio  keeps  increasing  marginally 
(for  central  core)  until  a  shear  strain  of  30%,  when  the  material  ceases  to  increase 
in  volume.  However  in  the  range  of  shear  strains  0.5%  to  30%,  the  shear  stress 
ratio  on  vertical  planes  for  the  whole  volume  marginally  decreases.  The  behavior 
for  the  central  core  and  whole  volume  are  contrary  to  each  other  and  this  could 
be  attributed  to  the  fact  that  there  is  more  mobility  of  the  particles  in  the  central 
core  than  at  the  boundaries.  Beyond  a  shear  strain  of  30%,  the  shear  stress  ratio 
on  vertical  planes  for  the  central  core  as  well  as  the  whole  volume  decreases 
sharply. 

(4)  The  shear  stress  ratio  on  horizontal  planes  reaches  a  peak  value  of  0.50  and 
decreases  at  larger  strains. 

The  observations  described  above  are  a  result  of  the  continuously  changing 
situation  concerning  the  stresses  (Figs.  4.22a,b).  While  a  constant  vertical  load  is 
applied,  the  normal  stresses  and  shear  stress  increases  and  reaches  a  peak  value 
at  a  shear  strain  of  0.04%,  when  a  mechanism  is  developed  (Fig.  4.22a).  At  large 
shear  strains,  the  normal  stress  on  horizontal  planes  (o^)  increases  marginally 
while  the  normal  stress  on  vertical  planes  (o^)  and  the  shear  stress  decreases. 
The  increase  or  decrease  in  the  values  of  the  stresses  in  the  central  core  at  the 
shear  strain  of  0.04%  relative  to  the  initial  vertical  load  conditions  is  as  follows. 
The  normal  stress  o^y  increases  by  11%  whereas  and  T^y  decreases  by  65% 


Shear  stress  (kPa)  Shear  stress  (kPa) 
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and  71%  respectively.  The  contribution  to  0^^^^  in  the  Hill's  averaging  principle 
comes  from  the  lateral  loads  on  the  vertical  boundaries.  When  the  shear  strain 
is  increased,  the  component  of  the  lateral  load  in  tire  direction  of  X-axis  decreases 
and  therefore  decreases.  Shear  stress,  decreases  due  to  the  decrease  in 
the  shear  force  as  it  can  be  seen  from  the  softening  of  the  shear  stress  ratio  curve 
for  the  bottom  boundary  in  Fig.  4.20b. 

The  stresses  computed  for  the  central  core  are  higher  than  that  for  the 
whole  volume.  The  values  are  compared  at  the  shear  strain  of  0.04%  relative  to 
the  initial  condition;  is  above  by  40%,  is  above  by  28%  ando^  is  above 
by  29%.  It  is  to  be  noticed  that  the  maximum  difference  is  in  the  shear  stress 
while  the  differences  in  the  normal  stresses  are  almost  the  same.  The  difference 
could  be  attributed  to  the  non-uniformity  in  the  contact  force  flow  paths  at  the 
top  left  comer  zone  and  the  bottom  right  comer  zone  (Figs.  4.14c,d,e)  as 
compared  to  the  central  core  where  the  load  paths  are  oriented  diagonally  and 
imiformly. 

In  a  granular  material,  failure  is  initiated  when  there  is  a  continuous  plane 
of  slippage /sliding  within  the  material.  This  stage  is  reached  when  the  shear 
stress  ratio  at  contacts  on  a  plane  exceeds  the  coefficient  of  friction  of  the  material. 
Thus,  in  a  granular  material,  failure  is  initiated  in  a  zone  of  maximum  shear  stress 
ratio.  By  referring  to  the  stress  strain  behavior  of  the  dense  packing  in  Fig.  4.21a, 
it  can  be  seen  that  the  slope  of  the  shear  stress  ratio  curve  for  vertical  planes  is 
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very  steep  compared  to  the  curve  for  horizontal  planes.  Thus,  the  vertical  planes 
reach  the  failure  state  at  a  much  faster  rate  than  the  failure  on  horizontal  planes. 

De  Josselin  De  Jong  (1971)  proposed  that  failure  in  simple  shear  could 
occur  either  by  sliding  on  horizontal  planes  or  sliding  on  vertical  planes  along 
with  rotation.  But,  the  sample  will  choose  the  latter  mode  with  least  resistance, 
if  both  modes  of  failure  are  equally  possible  for  the  prevailing  boundary 
conditions.  The  results  from  the  LTM  concord  with  De  Josselin  De  Jong's 
hypothesis. 

The  initial  shear  modulus  for  the  dense  packing  can  be  evaluated  from  Fig. 
4.20a  in  which  the  elastic  limit  appears  to  have  reached  at  a  shear  strain  of  0.0002. 
The  shear  modulus  is  evaluated  as  287  Mpa.  Assuming  a  Poisson's  ratio  of  0.20, 
the  initial  tangent  elastic  modulus  is  evaluated  as  120  MPa.  This  value  concords 
with  the  range  of  the  typical  initial  elastic  modulus  for  dense  sands  (100-200  MPa) 
given  by  Bowles  (1988). 

4.8.6  Failure  mechanisms 

During  the  shearing  process,  the  load  paths  transform  in  such  a  way  that 
contact  forces  are  concentrated  in  stiff  chains  of  particles  (Figs.  4.14a  through 
4.14h).  The  shear  stress  ratios  at  contacts  along  these  chains  are  far  less  than  'p' 
and  thus  these  contacts  never  slide.  Sliding  contacts  can  be  seen  only  in  relatively 
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unloaded  directions  or  between  the  stiff  chains  of  particles.  From  the  Hill's 
averaging  principle,  the  angle  made  by  the  major  principal  stress  plane  'i|;'  is 
computed  as  30®,  at  a  shear  strain  of  0.04%,  when  a  mechanism  is  reached. 
Based  on  this  angle  the  directions  of  the  major  and  minor  principal  stresses  are 
shown  in  Fig.  4.10.  It  was  earlier  explained  that  the  non-sliding  contacts  forming 
stiff  chains  of  particles  are  along  the  line  of  contacts  parallel  to  line  AA  (Figs. 
4.14d-h)  and  that  the  line  of  contacts  along  horizontal  rows  of  contacts  along  C2- 
C2,  D2-D2  etc.  and  the  line  of  contacts  parallel  to  line  BB  become  bondless.  In 
Fig.  4.10,  it  is  shown  that  the  direction  of  minor  principal  stress  makes  an  angle 
of  30®  with  the  above  mentioned  line  of  contacts  where  debonding  takes  place. 
Therefore,  it  can  be  concluded  that  as  a  material  deforms,  the  material  chooses  a 
mode  of  failure  in  such  a  way  that  non-sliding  contacts  with  strong  bonds  are 
developed  in  the  major  principal  stress  direction  and  the  material  chooses  the 
minor  principal  stress  direction  as  the  preferential  direction  for  losing  bonded 
contacts.  The  contacts  between  the  major  and  minor  principal  stress  directions 
remain  as  sliding  contacts  or  as  "stick  and  slip"  contacts. 

Cundall  and  Strack  (1982  and  1983)  performed  numerical  tests  on 
assemblies  of  particles  using  DEM.  The  above  mentioned  observations  from  the 
lattice  type  model  were  observed  by  them  in  the  results  from  DEM.  Oda  and 
Konishi  (1974)  performed  simple  shear  tests  on  an  assembly  of  photoelastic  disks. 
The  frequency  distributions  of  the  contact  normals  are  given  in  Fig.  4  of  their 
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paper.  The  photoelastic  experiments  also  revealed  that  as  the  sample  deforms  in 
simple  shear,  the  percentage  of  contacts  in  the  major  principal  stress  direction 
increases  and  the  percentage  of  contacts  in  the  minor  principal  stress  direction 
decreases.  Thus  the  observations  from  numerical  tests  using  DEM  and  the 
photoelastic  experiments  concord  with  the  observations  in  the  lattice  type  model. 
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CHAPTER  5 

CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  Conclusions 

In  this  research,  a  lattice  type  model  that  can  describe  the  micro  mechanical 
behavior  of  particulate  media  has  been  formulated  and  numerical  tests  have  been 
conducted  on  dense  and  loose  assemblies  of  disks.  The  particulate  media  is 
transformed  into  a  lattice  comprising  bars  linking  the  contacts  of  the  particles  and 
a  truss  analysis  is  then  conducted  with  appropriate  joint  constraints  to  determine 
the  response  of  the  particulate  media  imder  loads.  The  important  features  of  the 
model  are: 

•  A  conceptually  simple  model  that  uses  established  methods  of  analysis. 

•  The  particulate  media  itself  is  simulated  by  bars  of  representative  stiffness 
unlike  other  existing  methods,  such  as  the  Discrete  Element  Method,  where 
the  particles  are  replaced  by  springs  at  the  contacts.  This  is  especially 
important,  if  one  wishes  to  consider  phenomena  such  as  particle  (or  grain) 
crushing  and  initial  particle  imperfections. 

The  results  from  tire  numerical  tests  on  the  ideal  assemblies  of  disks  in  the 
loosest  and  densest  packing  established  the  robustness  and  potential  of  the  lattice 
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type  model  to  provide  insights  into  the  micro  mechanical  response  of  particulate 
media.  Among  the  characteristics  exhibited  by  the  assemblies  of  disks  are: 

•  Deformation  occurs  through  groups  or  clusters  of  particles  that  slide, 
rotate,  subdivide  and  regroup.  This  is  consistent  with  observations  from 
tests  on  real  granular  materials  and  images  captured  by  photoelastic  tests. 

•  Even  for  an  ideal  assembly  of  disks,  the  non-ideal  simple  shear  strain 
deformation  introduced  non-uniformities  in  stresses  that  result  in  cluster 
formations. 

•  Certain  established  patterns  of  behavior  in  particular  media  such  as 
compression  in  loose  material  and  dilation  in  dense  material  are  predicted 
by  the  lattice  t5^e  model. 

•  Using  Hill's  averaging  stress  principle,  the  discrete  set  of  contact  forces 
was  transformed  into  average  stresses  over  a  selected  region  of  the  dense 
assembly  of  disks  to  determine  constitutive  relationships.  The  results  of 
the  lattice  type  model  were  consistent  with  constitutive  relationships 
observed  in  real  particulate  assemblies. 

5.1  Recommendations 

The  lattice  type  model  is  still  in  its  infancy  and  needs  further  development 
to  assess  whether  it  can  be  a  viable,  main  stream  model.  The  immediate  needs 
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are: 

1.  Inclusion  of  additional  contacts  during  deformation.  This  is  necessary  to 
simulate  the  formation  of  additional  contacts  when  a  particle  moves  into  a  new 
position  and  meets  more  neighbors  than  in  its  original  position. 

2.  A  comparative  study  with  the  Discrete  Element  Method. 

3.  A  parametric  study  using  random  packing,  random  sizes,  irregular  shapes,  bar 
stiffness,  constraints  for  particle  crushing,  and  initial  particle  imperfections. 

4.  Application  of  the  model  to  real  particulate  systems  such  as  sands.  Here,  an 
image  of  the  intact  sand  can  be  taken  and  then  digitized.  The  data  can  then  be 
used  to  transform  the  sand  structure  into  a  lattice  from  which  the  analysis  imder 
the  desired  loading  can  be  conducted. 

5.  A  three-dimensional  version  of  the  model  should  be  attempted  since  the 
analysis  would  be  similar  to  a  space  truss. 

6.  The  lattice  type  model  can  be  extended  to  handle  dynamic  loading. 
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APPENDIX  - 1 


A  -  1  :  Displacements  in  an  elastic  disk 


Fig.  Al.l  shows  an  elastic  disk  acted  upon  by  an  equal  and  opposite  load 
P  along  a  chord  of  the  circular  cross-section  of  the  disk. 


P 


Fig.  Al.l  Elastic  disk  imder  load  along  a  chord 
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The  expressions  for  stresses  in  this  elastic  disk  are  given  in  Frocht  (1948).  These 
expressions  were  derived  from  the  original  solutions  of  Michell  (1900  and  1902) 
for  stresses  in  an  elastic  half  space.  The  normal  stresses  and  shear  stress  O^,  0^ 
and  T^y  shown  in  Fig.  3.5  are  given  by  Equations  (Al.l),  (A1.2)  and  (A1.3) 
respectively. 


-2P  {y^^y){^-^^? 

Q  =  -  [  -  +  -  -  C  ] 

nf  / 
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where  t  is  the  thickness  of  the  disk,  d  is  the  diameter  of  the  disk,  a  is  the  angle 

between  OA  and  AB,  C  =  =  (x-x.,)^  +  (y^-y)^  and 

d 

rl  =  {x-x^^  +  (y-y^)^  =  (x-x.,)^  +  (y+y^)^.  From  the  above  expressions  for 
stresses,  the  displacement  gradients  can  be  obtained  from  the  expressions 


du 

dx 


I  (O,  -  V  op 


dv 

dy 


=  =  i  -  V  a,) 


(A1.4) 


(A1.5) 
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where  u  and  v  are  the  x  and  y  components  of  the  displacements,  are  the 
engineering  strains  in  the  x  and  y  directions  and  E  and  v  are  the  elastic  modulus 
and  Poisson's  ratio  of  the  material  of  the  disk. 

The  displacements  u  and  v  at  a  general  point  (xp,  5rp)  are  obtained  by  the 
integration  of  the  displacement  gradients  between  this  point  and  a  point  where 
u  and  V  displacements  are  zero.  The  displacement  v  is  zero  at  all  the  points 
along  x-axis.  However,  a  point  where  the  displacement  u  is  zero  cannot  be 
determined  in  this  disk  problem  as  there  is  a  rigid  body  motion.  Therefore,  only 
relative  displacement  of  the  point  (xp,  yp)  witii  respect  to  the  centroid  of  the  disk 
can  be  determined.  The  equation  of  the  line  from  origin  (centroid)  to  the  point 
(xp,  yp)  is  obtained  as  y  =  mx,  where  m  =  yp/xp.  The  displacements  u  and  v  are 
given  by 


u  = 


xp 

f 


dU 

dx 


I  dx 

y  =  mx 


V  = 


yp 

! 


dv 

dy 


(A1.6) 


(A1.7) 


The  displacement  of  point  A,  with  respect  to  the  centroid  O  is  obtained  by 
integrating  the  above  expressions  from  (0,  0)  to  (x.,,  y^).  However,  this  will  lead 
to  singularity  and  thus,  theoretically,  it  is  not  possible  to  determine  the 
displacement  of  A  relative  to  the  centroid  O.  Therefore,  the  displacement  at  A  is 
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evaluated  as  a  sum  of  the  displacement  of  point  D  where  the  material  yields 
along  the  line  OA  and  the  product  of  the  displacement  gradient  at  point  D  and 
the  length  of  DA.  It  can  be  seen  from  Equations.  (Al.l)  and  (A1.2)  that  the 
stresses  are  singular  at  point  A.  Therefore,  a  point  D  at  (Xg,  /g),  where  the 
material  yields,  is  obtained  by  setting  the  expression  for  the  major  principal 
stress  at  point  D  equal  to  the  yield  stress  of  the  material.  The  major  principal 
stress  is  given  by 


0+0 

o,  =  + 


'  V  +  T 


xy 


(A1.8) 


Since,  O^,  O^  and  are  functions  of  x  and  y,  the  stresses  at  point  D  are 
obtained  by  substituting  (Xg,  /g)  for  (x,  y)  in  the  expressions  for  the  stresses. 
Therefore,  the  expression  for  O.,  which  is  set  equal  to  the  yield  stress,  will  be  a 
function  of  (Xg,  /g) .  Since  (Xg,  yg)  is  a  point  along  OA,  substituting  for  =  m 
(where  m  =  slope  of  OA  =  y^/x^ )  will  leave  the  expression  for  o.,  in  terms  of  Xg 
and  thus  Xg  is  evaluated.  Then,  /g  is  evaluated  from  /g  =  /77  Xg.  Knowing 
(Xg,  yg),  the  expressions  for  evaluating  the  displacements  at  point  A  are  given  by 


o,  =  Uo  +  I  (X,  -  x,) 

(X3.y3) 


(A1.9) 


dv 


-  Yu  \ 

(X3.y3) 


(Al.lO) 
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where  u^,  and  are  the  x  and  y  components  of  the  displacements  at  A 

and  D  respectively.  Uq,  are  evaluated  from  Equations.  (A1.6)  and  (A1.7) 
respectively. 


165 


REFERENCES 


Adley,  M.  D.  and  Sadd,  M.  H.  (1992).  "Continuum  models  for  materials  with 
lattice-like  microstructure",  in  Computers  &  Structures  Vol.  43,  No.  1,  pp.  13-18. 

Bardet,  J.  P.  And  Proubet.  J.  (1991a).  "An  adaptative  relaxation  technique  for  the 
statics  of  granular  materials".  Computers  and  Structures,  39  (3/4);  221-229. 

Bardet,  J.  P.  and  Proubet,  J.  (1991b).  "A  numerical  investigation  of  the  structure 
of  persistent  shear  bands  in  granular  media".  Geotechnique,  41(4),  pp.  599-613. 

Bardet,  J.  P.  and  Proubet,  J.  (1992).  "A  shear  band  analysis  in  idealized  granular 
materials".  Journal  of  Engineering  mechanics,  ASCE,  118(2),  pp.  397-415. 

Bathurst,  R.  J.  and  Rothenburg,  L.  (1988).  "Micro-mechanical  aspects  of  isotropic 
granular  assemblies  with  linear  contact  interactions",  in  J.  Appl.  Mech,  Vol.  55,  pp. 
17-23. 


Bishop,  A.W.  (1973).  "Shear  strength  parameters  for  imdisturbed  and  remoulded 
soil  specimens".  In  Stress  Strain  Behavior  of  Soils,  Edited  by  R.H.G.  Parry,  pp.  3- 
58. 

Boardman,  W.G.  (1990).  "Wave  propagation  in  granular  media  simulated  by 
elastic  networks".  Master's  Thesis,  Dept.  Of  Mechanical  Engineering  and  Applied 
Mechanics,  Univ.  Of  Rhode  Island,  U.S.A. 

Bowles,  J.E.  (1988).  Foundation  analysis  and  design.  Me  Graw  Hill. 

Burman,  B.C.  (1974).  "Development  of  a  numerical  model  for  discontinua",  Aust. 
Geomech.J.G4,  No.  1,  pp.  13-22. 

Dragoo,  B.  (1995).  "Two-dimensional  photoelastic  modelling  of  granular  media 
using  image  processing",  M.  S.  thesis.  Department  of  Civil  Engineering  and 
Engineering  mechanics.  University  of  Arizona,  Tucson,  AZ. 

Budhu,  M.  (1979).  "Simple  shear  deformation  of  sands",  Ph.  D.  thesis,  Cambridge 
University,  Cambridge,  United  Kingdom. 

Budhu,  M.  (1988).  "Failure  state  of  a  sand  in  simple  shear",  in  Can.  Geotech.  J.,  25, 
pp.  395-400. 


166 


Budhu,  M.,  Ramakrishnan,  S.  and  Frantziskonis,  G.  (1994).  "Mechanics  of 
particulate  media  -  A  lattice  type  approach".  Proceedings;  Workshop  on 
Mechanics  and  Statistical  Physics  of  Particulate  Materials,  Institute  for  Mechanics 
and  Materials,  University  of  California,  San  Diego,  CA,  June  8-10, 1994. 

Chang,  C.S.  and  Misra,  A.  (1989a).  "Theoretical  and  experimental  study  of  regular 
packings  of  granules".  Journal  of  Engineering  Mechanics,  Vol.  115,  No.4,  pp.  704- 
720. 

Chang,  C.  S.  and  Misra,  A.  (1989).  "Computer  simulation  and  modelling  of 
mechanical  properties  of  particulates".  Computers  and  Geotechnics,  4,  pp.  269-287. 

Chang,  C.  S.  and  Liao,  C.  L.  (1990).  "Constitutive  relation  for  a  particulate 
medium  with  the  effect  of  particle  rotation",  Int.  J.  Solids  and  Structures,  Vol.  26, 
No.  4,  pp.  437-453. 

Chang,  C.  S.  (1990).  "Strain  tensor  and  deformation  for  granular  material".  Journal 
of  Engineering  Mechanics,  ASCE,  Vol.  116,  No.  4,  790-804. 

Chang,  C.  S.  and  Misra,  A.  (1990).  "Packing  structure  and  mechanical  properties 
of  granulates".  Journal  of  Engineering  mechanics,  ASCE,  Vol.  116,  No.  5,  pp.  1077- 
1093. 

Cundall,  P.  A.  and  Strack,  O.  D.  L.  (1979a).  "A  discrete  numerical  model  for 
granular  assemblies,"  Geotechnique,  29,  pp.  47-65. 

Cimdall,  P.  A.  and  Strack,  O.  D.  L.  (1979b).  "The  distinct  element  method  as  a  tool 
for  research  in  granular  media".  Report  to  NSF  concerning  grant  ENG  76-20711, 
Part  2,  Dept.  Civ.  Min.  Engrg.  Univ.  Minnesota. 

Cundall,  P.  A.  and  Strack,  O.  D.  L.  (1979c).  "The  development  of  constitutive  laws 
for  soil  using  the  distinct  element  method".  Numerical  Methods  in  Geomechanics, 
Aachen,  ed.  W.  Wittke,  pp.  289-298. 

Cundall,  P.  A.,  Drescher,  A.,  and  Strack,  O.  D.  L.  (1982).  "Numerical  experiments 
on  granular  assemblies;  Measurements  and  observations".  Proceedings,  lUTAM 
Conference  on  Deformation  and  Failure  of  Granular  Materials,  Delft,  Edited  by 
P.  A.  Vermeer  and  H.  J.  Luger,  A.  A.  Balkema/  Rotterdam,  pp.  355-370. 

Cundall,  P.A.  and  Strack,  O.D.L.  (1983).  "Modelling  of  microscopic  mechanisms 
in  granular  material".  Mechanics  of  granular  materials  edited  by  J.  T.  Jenkins  and 


167 


M.  Satake,  Elsevier,  pp.  137-149. 

Cundall,  P.  A.  (1988).  "Computer  simulations  of  dense  sphere  assemblies", 
Micromechanics  of  granular  materials,  Ed.  by  M.  Satake  and  J.  T.  Jenkins,  Elsevier, 
Amsterdam,  pp.  113-123. 

Cundall,  P.A.  (1989).  "Numerical  experiments  on  localization  in  frictional 
materials",  Ingenieur-Archiv,  Vol.  59,  No.  2,  pp.  148-159. 

Cundall,  P.  A.,  Jenkins,  J.  T.  and  Ishibashi,  I.  (1989).  "Evolution  of  elastic  moduli 
in  a  deforming  granular  assembly".  Powders  and  Grains,  Ed.  by  J.  Biarez  and  R. 
Gourves,  Balkema,  pp.  319-322. 

Dai,  H.  and  Frantziskonis,  G.  (1994).  "Heterogeneity,  spatial  correlations,  size 
effects  and  dissipated  energy  in  brittle  materials".  Mechanics  of  materials, 
Elsevier,  pp.  103-118. 

Dantu,  P.  (1957).  "Contribution  a  I'etude  mecanique  et  geometrique  des  milieux 
pulverulents",  Proc.  4th  Int.  Conf.  Soil  Mech.  Foimd.  Engrg,  London:  144-148. 

De  Josselin  De  Jong,  G.  and  Verruijt,  A.  (1969).  "Etude  photo-elastique  d'un 
epilement  de  disques",  Cah.  Grpe  Fr.  Etud.  Rheol.  2  :73  -  86. 

De  Josselin  De  Jong,  G.  (1971).  Stress-strain  behavior  of  soils.  In  Proceedings, 
Roscoe  Memorial  Symposium,  Edited  by  R.H.G.  Parry,  U.K.  pp.  258-261. 

Dow,  J.  O.,  Su,  Z.  W.,  Feng,  C.  C.  and  Bodley,  C.  (1985).  "Equivalent  continuum 
representation  of  structures  composed  of  repeated  elements",  AIAA  Journal,  23, 
pp.  1564-1569. 

Deresiewicz,  H.  (1958).  "Mechanics  of  granular  material",  Advd.  appl.  Mech.  5,  pp. 
233-306. 

Drescher,  A.  and  De  Josselin  De  Jong,  G.  (1972).  "Photoelastic  verification  of  a 
mechanical  model  for  the  flow  of  a  granular  material",  J.  Mech.  Phys.  Solids,  Vol. 
20,  pp.  337-351. 

Duncan,  J.  M.  and  Chang,  C.  Y.  (1970).  "Nonlinear  analysis  of  stress  and  strain  in 
soils",  J.  Of  Soil  Mech.  and  Found.  Div.,  ASCE,  Vol.  96,  No.  SM5,  pp.  1629-1653. 

Feda,  J.  (1982).  "Mechanics  of  particulate  materials",  Elsevier  Pub  Company. 


168 


Frantziskonis,  G.,  Karpur,  P.,  Matikas,  T.  E.,  Krishnamurthy,  S.  and  Jero,  P.  D. 
(1994).  "Fiber-matrix  interface  -  information  from  experiments  via  simulation". 
Composite  structures,  Vol.  29,  pp.  231-247. 

Frocht,  M.  M.  (1948).  "Photoelasticity"  Vol.  2,  John  Wiley  &  Sons,  New  York. 

Gill,  J.  J.  (1993).  "The  microstructural  response  of  granular  soil  under  uniaxial 
strain".  Final  report  to  AFOSR,  Phillips  Laboratory,  Air  Force  Material  Command, 
Kirtland  AFB,  NM. 

Hansen,  A.,  Roux,  S.,  and  Herrmann,  H.  J.  (1989).  J.  Phys.  France  50,  733. 

Herrmann,  H.J.  and  Roux,  S.  (1990).  "Statistical  models  for  the  fracture  of 
disordered  media,  North-Holland. 

Hertz,  H.  (1882).  "Uber  die  Behnmgsfester  Elastischer  Korper",  J.  Reine  Angew. 
Math,  Vol.  92,  No.  7,  pp.  156-171. 

Hill,  R.  (1963).  "Elastic  properties  of  reinforced  solids:  Some  theoretical  principles". 
Journal  of  mechanics  and  physics  of  solids,  Vol.  11,  No.  5,  pp.  357-372. 

Hrennikoff,  A.  (1941).  "Solution  of  problems  of  elasticity  by  the  framework 
method".  Journal  of  Applied  Mechanics,  ASME,  A-169-175. 

Ishibashi,  I.  and  Chen,  Y.C.  (1988).  "Dynamic  shear  moduli  and  their  relationship 
to  fabric  of  granular  materials".  Micromechanics  of  granular  materials,  Ed.  by  M. 
Satake  and  J.  T.  Jenkms,  Elsevier,  Amsterdam,  pp.  95-102. 

Johnson,  K.L.  (1985).  Contact  mechanics,  Cambridge  University  Press. 

Kishino,  Y.  (1988).  "Disc  model  analysis  of  granular  media".  Micromechanics  of 
granular  materials,  Ed.  by  M.  Satake  and  J.  T.  Jenkms,  Elsevier,  Amsterdam,  pp. 
143-152. 

Kollar,  L.  and  Hegedus,  1.  (1985).  "Analysis  and  design  of  space  frames  by  the 
continuum  method",  Elsevier. 

Konishi,  J.  and  Naruse,  F.  (1988).  "A  note  on  fabric  in  terms  of  voids".  In 
Micromechanics  of  granular  materials,  Ed.  by  M.  Satake  and  J.  T.  Jenkins,  Elsevier, 
Amsterdam,  pp.  39-46. 


169 


Konishi,  Oda,  M.,  and  Nemat-Nasser,  S.  (1983),  "Induced  anisotropy  in 
assemblies  of  oval  cross-sectional  rods  in  biaxial  compression".  Mechanics  of 
granular  materials  edited  by  J.  T.  Jenkins  and  M.  Satake,  Elsevier,  pp.  31-39. 

Konishi,  J.  and  Naruse,  F.  (1988).  "A  note  on  fabric  in  terms  of  voids". 
Micromechanics  of  granular  materials,  Ed.  by  M.  Satake  and  J.  T.  Jenkins,  Elsevier, 
Amsterdam,  pp.  39-46. 

Krajcinovic,  D.  and  Silva,  M.A.  (1982).  "Statistical  aspects  of  the  continuous 
damage  theory",  Int.  J.  Solids  and  Structures,  Vol.  18,  No.  7,  pp.  551-562. 

Krajcinovic,  D.  and  Basista,  M.  (1991).  "Rupture  of  central-force  lattices  revisited", 
in  J.  Phys.  1 , 1,  pp.  241-245. 

Luerkens,  D.  W.  (1991).  "Theory  and  application  of  morphological  analysis",  CRC 
press. 

Mantell,  C.L.  (1958).  Engineering  materials  handbook,  McGraw  Hill. 

Mindlin,  R.D.  (1949).  "Compliance  of  Elastic  Bodies  in  Contact,"  J.  Appl.  Mech., 
Sept.,  259-268. 

Mindlin,  R.  D.  and  Deresiewicz,  H.  (1953).  "Elastic  spheres  in  contact  imder 
varying  oblique  forces",  J.  Appl.  Mech.,  20,  327-344. 

Michell,  J.  H.  (1900).  Proc.  London  Math.  Soc.,  Vol.  32,  p.  44. 

Michell,  J.  H.  (1902).  Proc.  London  Math.  Soc.,  Vol.  34,  p.  134. 

Misra,  A.  (1990).  "Constitutive  relationships  for  granular  solids  with  particle 
slidings  and  fabric  changes",  Ph.  D  dissertation.  University  of  Massachusetts, 
Amherst,  MA. 

Nemat-Nasser,  S.  and  Mehrabadi,  M.  (1983).  "Stress  and  fabric  in  granular 
masses",  in  Mechanics  of  Granular  materials  edited  by  Jenkins,  J.T  and  Satake,  M., 
Elsevier,  pp.  1-8. 

Nemat-Nasser,  S.  and  Mehrabadi,  M.  (1984).  "Micromechanically  based  rate 
constitutive  descriptions  for  granular  materials",  in  Mechanics  of  Engineering 
Materials,  Edited  by  C.S.Desai  and  R.H.Gallagher,  John  Wiley  &  Sons  Ltd. 


170 


Newmark,  N.M.  (1949),  "Numerical  meti\ods  of  analysis  in  bars,  plates  and  elastic 
bodies,"  in  'Engineering'  edited  by  L.  E.  Grinter,  Macmillan. 

Noor,  A.  K.  and  Russell,  W.  C.  (1986).  "Anisotropic  continuum  models  for 
beamlike  lattice  trusses".  Comp.  Meth.  Appl.  Mech.  Engrg.,  57,  pp.  257-277. 

Oda,  M.  (1972a).  "Initial  fabrics  and  their  relations  to  mechanical  properties  of 
granular  material".  Soil  and  Foundations,  Journal  of  the  JSSMFE,  Vol.  12,  No.  1, 
pp.  17-36. 

Oda,  M.  (1972b).  "The  mechanism  of  fabric  changes  during  compressional 
deformation  of  sand".  Soils  and  Foundations,  Journal  of  the  JSSMFE,  Vol.  12,  No. 

2,  pp.  1-18. 

Oda,  M.  (1972c).  "Deformation  mechanism  of  sand  in  triaxial  compression  tests". 
Soils  and  Foimdations,  Journal  of  the  JSSMFE,  Vol.  12,  No.  4,  pp.  46-63. 

Oda,  M.  and  Konishi,  J.,  (1974),  "Microscopic  Deformation  of  Granular  Material 
in  Simple  Shear,"  Soils  and  Foimdations,  14:  25-38. 

Oda,  M.,  (1978),  "Significance  of  fabric  in  granular  mechanics".  Proceedings,  U.S.- 
Japan  seminar  on  Continuum-Mechanical  and  Statistical  Approaches  in  the 
Mechanics  of  Granular  Materials,  Ed.  By  S.C.  Cowin  and  M.  Satake,  Gakujutsu 
Bunken  Fukyu-kai,  pp.  7-26. 

Oda,  M.,  Nemat-Nasser,  S.  and  Mehrabadi,  M.M.  (1982),  "A  statistical  study  of 
fabric  in  a  random  assembly  of  spherical  granules",  Int.  J.  Anal.  Methods  in 
Geomech.  6,  pp.  77-94. 

Paikowsky,  S.G.  and  DiRocco,  K.J.  (1993).  "Image  analysis  of  interparticle  contact 
modeling".  Digital  Image  Processing:  Techniques  and  Applications  in  Civil 
Engineering:  Proceedings  of  a  Conference,  Kona,  Hawaii,  Feb.  28  -  Mar.  5. 

Parikh,  P.  V.  (1967).  "The  shearing  behavior  of  sand  under  axis5nnmetric  loading", 
Ph.  D  thesis.  University  of  Manchester. 

Ramakrishnan,  S.  (1996).  "Mechanics  of  Particulate  Media  -  A  Lattice  Type 
Approach",  Ph.D  dissertation.  Department  of  Civil  Engineering,  University  of 
Arizona,  Tucson,  U.S.A,  In  progress. 

Roscoe,  K.  H.,  (1953),  "An  apparatus  for  the  application  of  simple  shear  to  soil 


171 


samples,"  Proceedings,  3rd  Int  Conference  on  Soil  Mechanics  and  Foundation 
Engineering,  Zurich,  Vol.  1,  pp.  186-191. 

Rothensburg,  L.  and  Bathurst,  R.  J.  (1992).  "Micromechanical  features  of  granular 
assemblies  with  planar  elliptical  particles".  Geotechnique,  42,  79-95. 

Sadd,  M.  H.,  Qui,  L.,  Boardman,  W.  G.,  Shukla,  A.  (1992).  "Modelling  wave 
propagation  in  granular  media  using  elastic  networks",  Int.  J.  Rock  Mech.  Min. 
Sci.  &  Geomech.  Abstr.,  Vol.  29,  No.  2,  pp.  161-170. 

Serrano,  A. A.  and  Rodriguez-Ortiz,  J.M.  (1973)  "A  contribution  to  the  mechanics 
of  heterogeneous  granular  media".  Proceedings,  Symposium  on  Plasticity  and  Soil 
Mechanics,  Cambridge,  U.K. 

Strack,  O.  D.  L.  and  Cundall,  P.  A.,  (1978).  "The  Distinct  Element  Method  as  a  tool 
for  Research  in  Granular  Media,"  Report  to  NSF  concerning  grant  ENG  76-20711, 
Part  1,  Dept.  Civ.  Min.  Engrg.  Univ.  Minnesota. 

Shukla,  A.  and  Damania,  C.  (1987).  "Experimental  investigation  of  wave  velocity 
and  dynamic  contact  stresses  in  an  assembly  of  disks".  Experimental  mechanics, 
30,  80-87. 


Subhash,  G.,  Nemat-Nasser,  S.,  Mehrabadi,  M.  M.,  Shodja,  H.  M.  (1991). 
"Experimental  investigation  of  fabric-stress  relations  in  granular  materials". 
Mechanics  of  materials,  Elsevier,  pp.  87-106. 

Thornton,  C.  and  Barnes,  D.  J.,  (1986).  "Computer  Simulated  Deformation  of 
Compact  Granular  Assemblies,"  ACTA  Mechanica,  64,  45-61. 

Thornton,  C.  And  Sun,  G.  (1993).  "Axisymmetric  compression  of  3D  polydisperse 
systems  of  spheres".  Powders  and  Grains,  Ed.  by  C.  Thornton,  Balkema,  pp.  129- 
134. 

Timoshenko,  S.P.  and  Goodier,  J.N.  (1970).  Theory  of  elasticity,  McGraw  Hill. 

Trollope,  D.  H.  and  Burman,  B.  C.,  (1980),  "Physical  and  Numerical  Experiments 
with  Granular  Wedges.",  Geotechnique,  30,  137-157. 

Wakabayashi,  T.  (1957).  "Photoelastic  method  for  determination  of  stress  in 
powdered  mass",  Proc.  7th  Jap.  Nat.  Congr.  Appl.  Mech.  pp.  153-192. 


PART  II 


PHOTOELASTIC  MODELING  OF  GRANULAR  MEDIA  USING  IMAGE 

PROCESSING 


by 


Muniram  Budhu 
Brian  Dragoo 
S.  Ramakrishnan 


A  report 
to  the 

National  Science  Foundation 
and 

Air  Force  Office  of  Scientific  Research 


June,  1996 


1 


TABLE  OF  CONTENTS 

CHAPTER  PAGE 

LIST  OF  FIGURES . 3 

ABSTRACT . 5 

CHAPTER  1 . 6 

1.0  Introduction . 6 

1.1  Objective . 6 

1.2  Scope . 6 

CHAPTER  2 

LITERATURE  REVIEW . 8 

2.0  Introduction . 8 

2.1  Photoelasticity . 8 

2.2  Photoelastic  Analyses  of  Contact  Forces . 11 

2.3  Micro-mechanisms  in  the  Deformation  of  Granular  Media . 15 

2.3.1  Flow  of  a  Granular  Assembly . 15 

2.3.2  Investigation  of  Granular  Materials  in  Simple  Shear . 18 

2.4  Digital  Image  Processing . 19 

2.4.1  Image  Processing  in  Photoelastic  Analysis . 19 

CHAPTER  3 

SYSTEM  SETUP . 21 

3.0  Introduction . 21 

3.1  The  Polariscope . 21 

3.2  Photosensitive  Material . 25 

3.2.1  Calibration  of  Material  Fringe  Constant . 25 

3.2.2  Calibration  of  P/Nd  Constant . 29 

3.3  Disks  Arrangements  and  Loading  System . 30 

3.4  Image  Processing  System . 33 

CHAPTER  4 

DEVELOPMENT  OF  CONTACT  FORCE  ANALYSIS  METHOD . 34 

4.0  Introduction . 34 

4.1  De  Josselin  de  Jong  and  Verruijt  Method  of  Force  Analysis . 34 

4.1.1  Distributed  vs.  Concentrated  Loads . 36 

4.1.2  Semi-infinite  Plane  vs.  Curving  Plane . 40 

4.1.3  Superimposition  of  Fringes . 42 

4.2  Application  of  Force  Analysis  Method . 45 

4.2.1  Contact  Force  Magnitude . 46 

4.2.2  Contact  Force  Direction . 49 

4.3  Semi-automation . 50 

CHAPTER  5 . 51 

RESULTS  AND  DISCUSSION . 51 

5.1  Behavioral  Observations . 51 


2 


TABLE  OF  CONTENTS  (cont.) 

5.1.1  Interdisk  Frictional  Behavior . 51 

5.1.2  Volumetric  Strain . 52 

5.1.3  Disks  Rotations . 55 

5.1.4  Load  Distribution . 57 

5.1.5  Uniformity  of  Boimdary  Stress . 63 

5.2  Comparison  With  the  Lattice  Model . 65 

CHAPTER  6 . 70 

CONCLUSIONS  AND  RECOMMENDATIONS . 70 

6.1  Conclusions . 70 

6.2  Recommendations . 72 

REFERENCES . 73 


3 


LIST  OF  FIGURES 

FIGURE  PAGE 

Figure  2.1  Photoelastic  isochromatic  fringes . 10 

Figure  2.2.  Photoelastic  image  of  two  disks . 12 

Figure  2.3.  Lines  of  constant  principal  stress  difference . 15 

Figure  2.4.  Schematic  diagram  of  load /deformation  mechanism . 16 

Figure  2.5.  Photoelastic  image  of  an  assembly  of  disks . 17 

Figure  3.1.  System  setup . 21 

Figure  3.2.  General  arrangement  of  a  lens  polariscope . 22 

Figure  3.3.  Diffused-light  polariscope . 24 

Figure  3.4.  Comparison  of  photoelastic  images . 24 

Figure  3.5.  Photoelastic  disk  in  diametrical  compression . 26 

Figure  3.6.  Diametrical  compression  load  vs.  fringe  order . 28 

Figure  3.7.  Concentrated  load  on  a  semi-infinite  plane . 29 

Figure  3.8.  (a)  Loose  and  (b)  dense  packing  extremes . 30 

Figure  3.9.  Mechanism  for  loading  and  shearing  disk  assembly . 31 

Figure  3.10.  Close-up  of  simple  shear  box . 32 

Figure  4.1.  Comparison  of  theoretical  and  photoelastic  stress  patterns . 35 

Figure  4.2.  Close-up  photoelastic  images  of  disk  contacts . 37 

Figure  4.3.  Disk  in  compression . 39 

Figure  4.4.  Disk  m  diametrical  compression . 40 

Figure  4.5.  Disk  in  contact  with  half-plane . 41 

Figure  4.6.  (a)  Loading  configurations  for  1.2"  disk . 43 

Figure  4.7.  Full-field  image  of  a  dense  packing . 46 

Figure  4.8.  Close-up  image  of  a  0.75-in.  diameter  disk . 47 

Figure  4.9.  Magnification  of  upper  left  contact  point . 48 

Figure  4.10.  Centroid  of  fringe . .! . 49 

Figure  5.1.  Photoelastic  disk  arrangements . 52 

Figure  5.2.  Volumetric  strain  vs.  shear  strain . 54 

Figure  5.3.  Rotation  of  grains . 55 

Figure  5.4.  Rotation  of  grains . 56 

Figure  5.5.  Rotation  of  grains . 57 

Figure  5.6.  Loose  packing  before  shear  strain  application . 58 

Figure  5.7a.  Photoelastic  image  of  loose  assembly  under  simple  shear . 59 

Figure  5.7b.  A  second  photoelastic  image  of  loose  assembly . 60 

Figure  5.7c.  A  third  photoelastic  image  of  loose  assembly . 61 

Figure  5.8.  Bookstack  analogy  of  simple  shear  deformation . 62 

Figure  5.9.  Boundary  nonuniformities . 64 

Figure  5.10.  Photoelastic  image  of  dense  packing  with  no  shear  strain . 65 

Figure  5.11.  Photoelastic  image  of  dense  packing  with  shear  strain . 66 


4 


FIGURE  PAGE 

Figure  5.12a.  Contact  force  vector  plot  of  dense  assembly  at  zero  shear 
strain  (Fig  5.10)  created  by  using  the  modified  version  of  de  Josselin  de 

Jong  and  Verruijt's  method  of  finding  contact  force  vectors . 67 

Figure  5.12b.  Contact  force  vector  plot  at  zero  applied  shear  strain  (Fig 

5.10)  predicted  by  the  lattice  type  model . 67 

Figure  5.13a.  Contact  force  vector  plot . 69 

Figure  5.13b.  Contact  force  vector  plot . 69 


ABSTRACT 


A  simple  shear  apparatus  was  designed  and  built  to  study  force 
distribution  within  an  assembly  of  photoelastic  disks  to  compare  with 
predictions  from  the  lattice  t5^e  model  described  in  Part  1.  The  principles  and 
methods  of  photoelasticity  were  used  to  gain  insights  into  how  micro-level 
parameters  affect  the  macro-level  behavior  of  a  particulate  media. 
Photosensitive  disks  were  arranged  in  the  simple  shear  apparatus  in  dense 
and  loose  packings  to  simulate  granular  assemblies.  Images  of  the 
displacements  and  deformation  of  the  photosensitive  disk  assemblies  were 
captured  with  a  video  camera,  and  an  image  processing  software  was 
modified  to  semi-automate  the  photoelastic  analysis.  Particle  deformation 
mechanisms,  volumetric  strain  relationships,  and  boundary  stress  patterns 
observed  during  simple  shear  testing  were  shown  to  be  consistent  with 
existing  theories  and  experimental  observations.  Load  redistribution  under 
simple  shear  was  observed  to  be  extremely  sensitive  to  frictional  micro¬ 
variations  on  disk  surfaces.  The  force  distributions  obtained  from  the 
photoelastic  analyses  were  consistent  with  predictions  of  the  lattice  type 
model  and  lends  credence  to  this  model  as  a  suitable  approach  to  study  the 
micromechanical  behavior  of  particulate  systems. 
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CHAPTER  1 

1.0  Introduction 

In  Part  I  of  this  report,  the  theoretical  basis  for  the  lattice  type  model 
(LTM)  has  been  described.  Numerical  experiments  were  conducted  to 
establish  the  procedures  of  and  predictions  from  the  model.  The  applicability 
of  LTM  depends  on  whether  its  predictions  are  consistent  with  results  from 
experiments  or  observed  behavior  of  particulate  media.  A  key  feature  of  LTM 
is  the  prediction  of  the  force  distribution  within  a  simulated  particulate 
assembly.  One  way  of  establishing  the  goodness  of  the  prediction  is  to 
determine  the  force  distribution  in  an  ideal  system  in  a  test  setup  and  then 
compare  the  test  results  with  the  model  prediction  for  the  ideal  system.  A 
test  method  that  offers  promise  is  photoelastic  experimental  method  and  it 
was  selected  to  evaluate  the  LTM. 


1.1  Objective 

The  objective  of  this  study  is  to  observe  and  determine  the  force 
distribution  and  the  behavior  of  packings  of  photoelastic  disks  in  simple 
shear  to  compare  with  predictions  from  the  lattice  type  model. 

1.2  Scope 

The  following  are  included  in  the  scope  of  this  work: 

•  Design  and  construction  of  a  simple  shear  box  capable  of  (a) 
applying  a  constant  vertical  load  to  a  two-dimensional  assembly  of 
disks,  (b)  deforming  the  assembly  under  simple  shear,  and  (c) 
transmitting  light  for  analysis  by  photoelasticity. 
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•  Develop  a  method  of  contact  force  analysis  using  image  processing. 

•  Perform  simple  shear  experiments  on  both  loose  and  dense 
assemblies  of  disks. 

•  Make  observations  for  general  comparison  with  current  theories 
and  observations  on  granular  media  under  simple  shear  loading. 

•  Compare  force  distribution  with  predictions  from  the  lattice-type 
model. 
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CHAPTER  2 
LITERATURE  REVIEW 


2.0  Introduction 

This  chapter  provides  general  information  on  prior  research  in  areas  of 
interest  to  the  present  work.  Areas  covered  include  photoelasticity, 
experimental  investigations  of  microscopic  mechanisms  in  the  deformation 
of  granular  materials,  photoelastic  investigations  of  contact  forces,  elastic 
analyses  of  contact  forces,  and  digital  image  processing  and  its  applications  in 
photoelastic  analysis. 

2.2  Photoelasticity 

Stress  analysis  has  proven  valuable  in  the  field  of  engineering.  The 
ability  to  determine  how  stresses  are  distributed  in  a  machine  or  structure 
under  a  given  set  of  loads  is  crucial  to  its  design.  The  most  widely  used 
method  for  stress  analysis  is  the  theory  of  elasticity  -  a  mathematical  theory 
that,  while  accurate,  becomes  increasingly  involved  as  geometry  and  loading 
conditions  become  more  complex.  It  is  often  difficult  for  engineers  to  find 
solutions  to  practical  problems  due  to  complex  boundary  conditions.  Hence, 
experimental  techniques  are  often  called  upon  to  aid  in  problem  solving. 

Photoelasticity  is  an  experimental  technique  that  takes  advantage  of  the 
optical  properties  of  a  transparent  material  to  quantify  the  distribution  of 
stress  within  that  material.  Machine  parts,  tools,  structures,  and  other  load- 
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bearing  objects  can  be  analyzed  for  stress  concentration  and  other  failure- 
inducing  phenomena  when  modeled  and  subjected  to  a  photoelastic  analysis. 
Photoelasticity,  as  a  science,  first  appeared  around  the  turn  of  the  century,  and 
by  the  1930s  it  was  developed  into  a  feasible  technique  for  stress  analysis  by 
Coker  and  Filon  at  the  University  of  London  (Frocht,  1941). 

A  model  specimen  cut  from  a  transparent,  photosensitive  material 
reacts  to  stress  by  undergoing  a  change  in  index  of  light  refraction  that 
Maxwell,  in  1853,  observed  to  be  proportional  to  the  stresses  induced  on  the 
specimen.  This  change  can  be  quantified  by  subjecting  the  specimen  to  a  field 
of  polarized  light  in  an  apparatus  known  as  a  polariscope  (see  Chapter  3).  The 
light  that  is  passed  through  the  specimen  is  extinguished  in  banded  regions 
called  fringes  (Fig.  2.1)  that  are  related  to  either  the  magnitude  or  direction  of 
the  principal  stresses.  If  the  polariscope  is  further  modified  (with  the  addition 
of  optical  plates),  the  light  field  can  become  circularly  polarized.  This 
technique  allows  the  fringes  related  to  the  principal  stress  directions 
(isoclinics)  to  be  eliminated,  and  the  remaining  stress  magnitude-related 
fringes  (isochromatics)  can  be  used  to  calculate  the  magnitude  of  the  principal 
stress  difference  (oi  -  02)  at  any  point,  which  is  equal  to  twice  the  maximum 
shear  stress  (Xmax)  at  that  point. 

The  fringes  develop  as  ordered  bands  as  the  stress  in  a  specimen  is 
increased,  with  fringe  order  (N)  zero  developing  first,  then  fringe  order  1,  etc. 
The  fringes  are  contours  of  equal  principal  stress  difference,  which  is  linearly 
related  to  the  fringe  order  by  a  the  stress-optic  law. 
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(2.1) 
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Figure  2.1  Photoelastic  isochromatic  fringes  on  a  disk  under  four  equal  radial  loads  (from 
Frocht,  1941).  Fringes  are  lines  of  constant  principal  stress  difference  and  are  ordered  as 
integers,  with  fringe  order  zero  the  first  to  develop.  Fringe  order  zero  at  the  center  reveals  an 
isotropic  point. 

where  fo  is  the  fringe  constant  of  the  photosensitive  material,  and  h  is  the 
thickness  of  the  specimen  (Dally,  1965).  If  f^  is  calibrated,  and  h  is  measured, 
the  principal  stress  difference  (and  hence  tmax)  can  be  read  directly  from  a 
photoelastic  image  for  any  point  in  a  specimen  falling  on  an  integered  fringe, 
and  can  be  interpolated  for  those  falling  between  fringes.  This  technique 
revolutionized  the  field  of  stress  analysis  because  it  allowed  machine  parts  or 
structures  with  shapes  that  were  too  complex  to  be  readily  handled  by 
elasticity  to  be  accurately  analyzed  for  design  purposes. 
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The  finite  element  method  (FEM)  and  other  numerical  techniques 
used  in  stress  analysis  for  design  have  since  replaced  photoelasticity  in  many 
fields  due  to  the  ease  with  which  boundary  conditions  can  be  changed  in 
numerical  models.  In  contrast,  the  study  of  different  conditions  by 
constructing  and  testing  several  photoelastic  specimens  is  time  consuming 
and  requires  a  high  degree  of  skill. 

However,  in  the  emerging  field  of  non-continuum  modeling, 
including  the  current  study  of  the  problem  of  stresses  flowing  in  a  granular 
assembly,  photoelasticity  can  be  used  to  determine  the  feasibility  of  proposed 
numerical  methods. 

2.2  Photoelastic  Analyses  of  Contact  Forces 

De  Josselin  de  Jong  and  Verruijt  (1969)  developed  a  technique  to  apply 
the  principles  of  photoelasticity  to  the  specific  task  of  finding  the  magnitude 
and  direction  of  a  contact  force  between  two  disks  (Fig.  2.2). 


Figure  2.2.  Photoelastic  image  of  two  disks  representing  an  intergranular  contact  point  in  a 
particulate  media. 

In  a  contact  stress  problem  on  a  linear  elastic  material,  the  geometry  of 
the  pattern  of  principal  stress  difference  that  develops  (as  represented  by  the 
photoelastic  fringes)  is  known  to  be  proportional  to  the  concentrated  force  P 
such  that: 

(T,-cT2  =  P-f(jc,y)  (2.2) 

where  f(x,y)  is  a  function  that  expresses  the  shape  of  the  isochromatic  fringes. 

Since  the  region  of  interest  in  the  disk  contact  stress  problem  is  close  to 
the  edge  of  the  disk,  de  Josselin  de  Jong  assumed  that  f(x,y)  could  be 
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approximated  by  a  concentrated  force  acting  on  a  semi-infinite  plane  as  given 
by  Frocht  (1948): 


2 

nhd 


(2.3) 


where  d  is  the  diameter  of  the  circle  that  best  approximates  the  fringe  passing 
through  both  the  contact  point  and  the  point  (x,y),  the  center  of  the  circle 
being  on  the  line  of  action  of  the  contact  force. 

By  rearranging  Eqns  2.1,  2.2,  and  2.3,  it  can  be  shown  that 

p 

—  =  constant  (2.4) 

nd 

for  any  contact  force  through  a  given  material  with  a  given  thickness.  This 
implies  that  if  this  constant  were  calibrated  by  measuring  the  diameter  of  a 
known  fringe  order  produced  by  a  known  concentrated  force  on  a  semi¬ 
infinite  plane,  the  value  of  the  force  through  any  contact  could  be  determined 
by  simply  measuring  the  diameter  of  a  known  fringe  order,  and  calculating 
the  magnitude  of  P  from  the  constant.  And,  as  mentioned  above,  the  force  is 
known  to  act  along  a  line  connecting  the  contact  point  and  the  center  of  the 
circle  which  best  approximates  the  measured  fringe. 

This  method  greatly  simplifies  the  photoelastic  analysis  for  this  specific 
purpose.  However,  since  the  assumption  that  the  disk  is  a  semi-infinite 
plane  only  holds  true  near  the  contact  area,  care  must  be  taken  in  the  analysis 
so  that  this  assumption  is  still  valid.  Care  must  also  be  exercised  in  using  the 
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de  Josselin  de  Jong  and  Verruijt  method  of  contact  force  determination  when 
measuring  fringes  very  near  the  contact  area. 

Poritsky  (1950)  developed  equations  for  stresses  and  deflections  of  two 
cylindrical  bodies  in  contact  under  both  normal  and  tangential  distributed 
loads.  These  equations,  which  are  expressed  in  terms  of  the  Airy  functions, 
can  also  be  plotted  as  curves  of  constant  principal  stress  difference  for  direct 
comparison  with  photoelastic  fringes.  By  examining  a  plot  of  Poritsky's 
constant  principal  stress  difference  curves  of  a  distributed  contact  force  on  a 
semi-infinite  plane,  it  can  be  seen  that  fringe  distortion  is  severe  in  the  region 
near  the  contact.  (Fig  2.3).  Frocht  showed  from  photoelastic  studies  that 
fringes  in  this  region  will  be  circles  for  a  concentrated  load. 

Since  de  Josselin  de  Jong  and  Verruijt  assume  in  their  analysis  that  the 
contact  load  is  concentrated  rather  than  distributed,  some  discussion  of  how 
to  measure  fringes  near  the  area  of  contact  is  necessary.  Detailed  discussion  of 
a  modified  technique  for  contact  force  and  magnitude  determination  based  on 
the  de  Josselin  de  Jong  and  Verruijt  analysis  is  presented  in  Chapter  4. 


15 


Figure  2.3.  Lines  of  constant  principal  stress  difference  for  an  oblique,  distributed  load 
on  a  semi-infinite  plane  (from  Poritsky,  1950).  The  circled  region  shows  pronounced 
fringe  deformation  due  to  the  distributed  load. 

2.3  Micro-mechanisms  in  the  Deformation  of  Granular  Media 
2.3.1  Flow  of  a  Granular  Assembly 

Drescher  and  de  Josselin  de  Jong  (1972)  studied  an  assembly  of 
photosensitive  particles  to  gain  insight  into  how  a  granular  media  flows 
under  load.  Their  experimental  setup  allowed  them  to  deform  the  assembly 
similar  to  a  simple  shear  deformation  (Fig  2.4). 
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Figure  2.4.  Schematic  diagram  of  load /deformation  mechanism 
used  by  Drescher  and  de  Josselin  de  Jong  (1972)  for  photoelastic 
investigation  of  stress  flow  through  a  granular  assembly. 


It  is  apparent  from  the  photoelastic  fringe  patterns  of  Drescher  and  de 
Josselin  de  Jong's  assembly  (Fig  2.5)  that  forces  are  transmitted  not  uniformly 
throughout  the  assembly,  but  by  chains  of  aligned  contact  forces.  These 
chains  then  form  relatively  rigid  blocks  of  particles  that  slide  relative  to  one 
another.  They  concluded  that  since  the  probability  of  developing  a  chain  in  a 
specific  location  is  small,  that  the  distances  between  chains  must  be  large;  at 
least  several  disk  diameters.  The  result  of  this  is  a  discrete,  rather  than 
continuous  force  distribution,  completely  contrary  to  the  continuum 
approach  commonly  taken  for  stress  determination  in  soil  mechanics 
problems. 
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Figure  2.5.  Photoelastic  image  of  an  assembly  of  disks,  randomly  arranged  (from  Drescher 
and  de  Josselin  de  Jong  (1972).  Disks  with  dark  areas  exhibit  compression  stresses  and  thus 
belong  to  the  instantaneous  structural  skeleton  of  the  assembly. 


A  similar  problem  exists  for  deformation.  When  just  one  contact 
yields  to  sliding,  which  occurs  when  the  angle  between  the  force  between  two 
disks  and  the  normal  to  the  contact  exceeds  the  allowable  friction  angle,  many 
contact  points  can  be  reoriented.  This  leads  to  a  drastic,  discrete  redistribution 
of  the  contact  forces.  Thus  the  movements  of  the  particles  are  discrete  in 
character,  which  is  again  contrary  to  the  approach  taken  for  deformations  in 
soil  mechanics. 

Of  course,  it  is  impractical  to  attempt  a  discrete  stress  and  strain  analysis 
when  making  soil  mechanics  computations.  The  authors  introduce  an 
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average  stress  and  strain-rate  tensor  evaluation  to  interpret  behavior  in  a 
granular  assembly.  Experimental  results  can  then  be  compared  with 
theoretically  predicted  behavior  based  on  a  mechanical  model.  The  main 
observations  of  Drescher  and  de  Josselin  de  Jong  are  that  the  division  of  the 
assembly  into  rigid,  sliding  blocks,  and  the  free  rotation  of  these  blocks,  are 
observed  in  practice. 

2.3.2  Investigation  of  Granular  Materials  in  Simple  Shear 

The  deformation  of  a  granular  media  was  simulated  by  Oda  and 
Konishi  (1974)  using  an  assembly  of  cylindrical  rods  made  from 
photosensitive  material  packed  at  random  into  a  simple  shear  apparatus. 
They  concluded  that  the  tendency  of  the  concentration  of  the  normal  vectors 
to  the  contact  planes  was  determined  primarily  by  the  mobilized  stress  ratio 
ilcSjx,  where  t  is  the  shear  stress  and  an  is  the  normal  stress.  By  plotting 
frequency  distributions  of  contact  normals,  they  noticed  that  the  preferred 
direction  of  these  normals  gradually  rotates  with  the  rotation  of  principal 
stress  axes  during  shear  deformation.  They  also  concluded  that  sliding  is 
confined  to  only  a  preferred  set  of  contacts  at  any  one  time,  implying,  like 
Drescher  and  de  Josselin  de  Jong,  that  particles  group  themselves  into  a 
system  of  rigid  block  that  slide  against  one  another. 
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2.4  Digital  Image  Processing 

Digital  imaging  is  a  widely  used  computational  technique  for  digitally 
representing  an  image  by  discretizing  it  and  assigning  numerical  values  to  its 
components.  Greyscale  images,  the  type  used  in  this  project,  are  broken  down 
into  picture  elements,  or  pixels,  that  have  three  dimensions:  an  x-coordinate, 
a  y-coordinate,  and  a  brightness  value.  The  brightness  value  is  a  number 
between  0  (black)  and  255  (white)  that  represents  a  relative  luminosity  of  the 
pixel  on  an  eight  bit  (2®)  scale.  Once  an  image  has  been  digitized,  the  three- 
dimensional  array  that  represents  it  can  be  readily  manipulated  by  imaging 
software.  Imaging  techniques  such  as  smoothing,  sharpening,  and  finding 
edges,  which  are  valuable  for  extracting  necessary  information  from  the 
image,  are  based  upon  manipulation  of  the  brightness  values  exhibited  by 
each  pixel  in  the  image. 

2.4.1  Image  Processing  in  Photoelastic  Analysis 

A  large  body  of  literature  has  been  produced  in  the  past  fifteen  years  on 
the  subject  of  automating  photoelastic  data  collection  using  digital  image 
processing.  These  include  Muller  and  Saackel  (1979),  Hu  et  al  (1983),  Umezaki 
et  al  (1988),  Patterson  (1988),  Gillies  (1988),  Chen  and  Taylor  (1989),  Ramesh  et 
al  (1991)  and  many  others.  However,  most  of  these  works  explain  complex 
processing  algorithms  for  such  techniques  as  fringe  separation,  fringe 
thirming,  etc.  Since  pre-packaged  imaging  software  (which  includes  these  and 
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other  processing  algorithms)  is  used  in  this  study,  it  is  unnecessary  to  detail 
these  works. 

Paikowsky  and  DiRocco  (1993)  have  also  developed  an  image  analysis 
system  exclusively  for  analyzing  photoelastic  images  of  two-dimensional 
modeling  of  interparticle  contacts.  They  cite,  as  advantages  to  using  image 
processing  to  analyze  these  images,  (1)  high  resolution  of  enlarged  portions  of 
images  and  (2)  the  possibility  of  high  frequency  image  capturing  for  motion 
analysis.  Disadvantages  include  high  capital  cost  and  time  consuming 
analysis. 

They  have  also  concluded  that  two  of  de  Josselin  de  Jong  and  Verruijt's 
most  critical  assumptions  lead  to  incorrect  measurement  of  contact  force 
magnitudes  from  photoelastic  images.  These  assumptions  are  as  follows:  (1) 
the  contact  area  of  the  particles  is  negligible  and  that  (2)  the  fringes  that  result 
from  contact  forces  are  circles  that  pass  through  the  contact  point  .  They  have 
thus  developed  their  own  technique  for  force  magnitude  and  orientation 
determination,  which  is  mentioned  but  not  detailed  in  Paikowsky  et  al  (1993). 
It  will  be  shown  (Chapter  4)  that  the  method  of  force  determination 
developed  by  de  Josselin  de  Jong  and  Verruijt  can  be  slightly  modified  to  give 
satisfactory  measurements  of  both  force  magnitudes  and  directions  using 
digital  image  processing  technology. 
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CHAPTER  3 
SYSTEM  SETUP 

3.0  Introduction 


Photoelastic  analyses  require  the  following:  a  polariscope,  a  load 
application  mechanism,  a  photosensitive  model,  and  system  for  acquiring 
and  processing  photoelastic  images  (Fig.  3.1). 


MONITOR 


(Fig.  3.2)  f  Macintosh  Quadra  840AV  1 

\NIH  Image  Processing  Software  J 

ugure  3.1.  System  setup  showing  polariscope  (see  Fig.  3.2  for  more  detail),  shear  box  with 
photosensitive  disks,  camcorder  for  image  capture,  computer  for  image  processing,  and  monitor 
for  image  analysis  and  display. 


3.1  The  Polariscope 

The  polariscope  (Fig.  3.2)  is  a  system  of  lenses  and  plates  aligned  to 
produce  a  field  of  circularly  polarized  light.  It  generally  consists  of  a  light 
source,  a  polarizer,  two  quarter-wave  plates,  and  an  analyzer.  The  light 
source  produces  white  light,  which  is  often  immediately  broken  down  into 
monochromatic  light  with  a  filter. 
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Figure  3.2.  General  arrangement  of  a  lens  polariscope. 


The  polarizer  polarizes  the  light  by  restricting  the  light  vector  to  single 
plane  vibration.  The  light  then  passes  through  the  first  quarter-wave  plate 
which  circularly  polarizes  the  light.  A  light  field  is  circularly  polarized  when 
the  tip  of  the  light  vector  describes  a  circular  helix  as  the  light  propagates. 
When  this  circularly  polarized  beam  passed  through  the  photosensitive 
model  (set  between  the  quarter-wave  plates),  which  behaves  as  a  doubly 
refracting  crystal,  it  is  resolved  into  two  equally  plane-polarized  components 
in  mutually  perpendicular  planes  determined  by  the  state  of  stress  in  the 
model.  The  second  quarter-wave  plate  de-circularizes  the  light,  giving  it 
plane  polarization  once  again. 

The  light  finally  reaches  the  analyzer,  where  it  is  either  fully 
extinguished  (dark  field)  or  passed  on  (light  field),  depending  on  the  desired 
effect.  All  photoelastic  tests  for  the  present  study  were  conducted  in  light 
fields. 

Additional  optical  devices  are  sometimes  employed  depending  on  the 
light  source  and  desired  final  image.  For  example,  partial  optical  mirrors  can 
be  used  in  such  techniques  as  fringe  sharpening  or  fringe  thinning  if 
additional  stress  information  is  required  from  a  specimen. 
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Two  polariscopes  were  used  for  this  work:  a  lens  polariscope  with  a 
small  viewing  area  and  a  diffused-light  polariscope  with  a  large  viewing  area. 
The  lens  polariscope  (the  type  shown  in  Fig.  3.2),  manufactured  by  Polarizing 
Instrument  Company,  is  used  for  the  capture  of  images  where  sharp  detail  is 
required.  This  polariscope  generally  requires  more  exact  lens  alignment,  plate 
rotation  position,  and  camera  focal  length,  which  can  make  experiments 
relatively  time  consuming. 

The  light  from  the  source  of  a  lens  polariscope,  emitted  through  a 
small  aperture,  is  subsequently  diffused  by  a  field  lens  before  entering  the 
polarizer,  and  is  re-collected  by  another  field  lens  before  entering  the  camera. 
The  viewing  area  of  this  polariscope  was  not  large  enough  to  capture  full- 
field  images  of  some  of  the  models  required,  so  a  second  polariscope  was 
needed. 

The  diffused-light  polariscope  (Fig.  3.3),  constructed  in  1975  by 
undergraduate  students  of  photoelasticity  at  the  University  of  Arizona,  has 
sufficient  viewing  area  for  all  models  used,  and  adjustments  generally  do  not 
require  the  same  degree  of  precision  as  the  lens  polariscope.  Because  this 
polariscope  uses  a  diffused  light  source,  it  tends  to  reduce  the  sharpness  and 
contrast  of  photoelastic  images.  Dark  fringes  appear  gray,  rather  than  black  as 
in  images  from  the  lens  polariscope  (Fig.  3.4).  Images  from  the  diffused-light 
polariscope,  because  of  their  inferior  contrast,  are  not  sufficient  for  the  high- 
resolution  analysis  required  to  collect  quantitative  data.  The  increased 
viewing  area  of  the  diffused-light  polariscope  is,  however,  taken  advantage  of 
for  qualitative  analysis. 
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Figure  3.3.  Diffused-light  polariscope. 


Figure  3.4.  Comparison  of  photoelastic  images  captured  from  lens  polariscope  (left)  and 
diffused  light  polariscope.  Higher  quality  fringe  contrast  is  available  from  the  lens 
polariscope,  which  is  used  to  capture  images  with  better  resolution  for  qualitative 
analysis. 
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3.2  Photosensitive  Material 


The  models  used  to  simulate  grains  in  a  particulate  assembly  were  0.75- 
in.  diameter  disks  machined  from  PSM-1,  a  photosensitive  polycarbonate 
material  manufactured  by  Vishay  Measurements  Group.  Vishay  advertises 
PSM-1  as  a  material  that  is  excellent  for  two-dimensional  model  work,  is  free 
of  creep  and  edge  effects,  has  excellent  photosensitivity,  is  easy  to  machine,  is 
non-brittle,  and  shows  excellent  transparency.  It  exhibits  the  following 
material  properties: 


Table  1.  PSM-1  Material  Properties 


Elastic  Modulus  E,  psi  t 

3.60  X 10^ 

Poisson's  Ratio  v  + 

0.38 

Material  Fringe  Constant  psi-in  J 

92.0 

P/Nd  Constant,  Ib/in-fringe 

148.0 

Thickness  h,  in. 

0.25 

t  Property  given  by  manufacturer 
^  Property  calibrated  in  laboratory 


Models  were  machined  with  an  end  mill  at  high  rpm  to  avoid 
overheating  of  model  boundaries,  which  could  produce  unwanted  residual 
stress  effects. 


3.2.1  Calibration  of  Material  Fringe  Constant 

To  calibrate  the  material  fringe  constant  fa  of  the  PSM-1  photosensitive 
material  used  in  these  tests,  the  method  suggested  by  Dally  (1965)  was 
employed. 
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It  is  known  from  the  mechanics  of  materials  that  the  principal  stress 
difference  along  the  horizontal  diameter  of  a  disk  loaded  in  diametrical 
compression  (Fig.  3.5)  is  given  by  the  following  equation: 


CTj  (J2 


SP  Z)^-4Z)V 
7thD(^D^+4x^f 


where  D  is  the  disk  diameter,  x  is  the  distance  along  the  horizontal  diameter 
measured  from  the  center  of  the  disk,  and  h  is  the  disk  thickness. 


Figure  3.5.  Photoelastic  disk  in  diametrical  compression 


It  is  also  known  from  the  stress-optic  law  that 
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where  N  is  the  fringe  order.  Setting  these  equations  equal  and  solving  for  the 
material  fringe  constant  f^  yields 


8P  £)'*-4DV 
nDN[D^+Ax^)'' 


(3.3) 


If  only  the  center  of  the  disk  is  considered,  x  =  0,  and  Eq.  3.3  reduces  to 


/.= 


8P 

ttDN 


(3.4) 


Note  that  the  value  of  f^  is  independent  of  the  thickness  of  the  disk. .  A  test 
was  arranged  where  the  load  on  the  disk  in  diametrical  compression  was 
applied  in  increasing  increments  while  the  fringe  order  at  the  center  of  the 
disk  was  read  for  each  load.  A  plot  of  P  versus  N  was  made,  and  an  average 
P/N  value  was  foimd  and  substituted  into  Eq.  3.4  to  find  a  value  of  fa- 

The  results  of  the  above  procedure  are  shown  in  Figure  3.6.  The 
material  fringe  constant  f^  for  PSM-1  was  found  to  be  92.0  psi-in.  To  find  the 
principal  stress  difference  (and  thus  the  maximum  shear  stress)  in  psi  at  any 
point  in  any  model  made  from  this  material,  this  value  of  f^  is  simply 
multiplied  by  the  fringe  order  N  represented  at  that  point  and  divided  by  the 
thickness  of  the  model  in  inches  (Eq.  3.2). 


Material  Fringe  Constant  Calibration:  P  vs.  N 


P/N  =  27.1 

Material  Fringe  Constant  =  92.0  psi-in 


0 

2 

4  6 

8 

10 

Fringe  Order,  N 

Figure  3.6.  Diametrical  compression  load  vs.  fringe  order  at  the 
center  of  the  disk  for  the  calibration  of  the  material  fringe  constant 
f<j.  Calibration  readings  were  taken  on  half-integer  fringe  intervals. 


29 


3.2.2  Calibration  of  P/Nd  Constant 

De  Josselin  de  Jong  and  Verruijt  (1969)  determined  that  a  concentrated 
load  P  acting  on  a  semi-infinite  plane  produces  circular  fringes  of  constant 
maximum  shear  stress  such  that  P/Nd  is  a  constant  for  any  given  material 
(where  N  is  the  order  of  a  circular  fringe  of  diameter  d).  This  constant  can  be 
easily  calibrated  for  the  PSM-1  material  by  capturing  an  image  of  a 
concentrated  load  on  a  semi-infinite  plane  (Fig  3.7)  and  measuring  the 
diameter  of  fringes  of  known  order. 


calibrated  from  this  image,  is  148  Ib/in-fringe  for  the  PSM-1  material. 


Once  this  constant  is  known  for  the  material,  the  magnitude  of  a  force 
through  a  contact  point  can  be  found  by  measuring  the  diameter  of  a  fringe  of 
known  order,  and  by  solving  for  P  from  the  constant.  A  detailed  explanation 
of  this  procedure  and  its  underlying  theory  is  given  in  Chapter  4. 
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3.3  Disks  Arrangements  and  Loading  System 

The  0.75-in.  disks  were  arranged  in  two  assemblies:  a  simple  cubic 
packing  to  simulate  a  loose  media  and  a  simple  stagger  to  simulate  a  dense 
media  (Fig.  3.8).  These  packings  illustrate  the  extremes  of  granular  packings; 
that  is,  the  disks  carmot  be  packed  any  looser  than  the  simple  cubic  packing 
nor  any  denser  than  the  simple  stagger. 


Figure  3.8.  (a)  Loose  and  (b)  dense  packing  extremes 


The  assemblies  were  placed  into  the  simple  shear  box  (Fig.  3.9)  which 
consists  of  a  base  that  moves  horizontally  by  turning  a  screw;  a  top  platen 
allowed  to  move  only  vertically,  and  two  side  walls  allowed  to  pivot  at  the 
top  (by  pin  connection)  and  bottom  (by  knife  edge)  to  retain  the  parallelogram 
shape  of  the  assembly  (Fig.  3.10).  The  assembly  was  loaded  vertically  with  a 
100-psi  pneumatic  ram  capable  of  applying  491  lb  to  the  top  platen,  which 
contacted  the  top  row  of  disks.  The  assembly  was  subject  to  simple  shear  by 
moving  the  base  horizontally.  The  top  platen  was  allowed  to  move  vertically 
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as  the  assembly  dilated  or  compressed  while  keeping  the  normal  load 
constant. 


Figure  3.9.  Mechanism  for  loading  and  shearing  disk  assembly. 


Shear  strains  were  calculated  by  measuring  horizontal  displacements  of  the 
base  with  a  dial  gauge,  and  volumetric  strains  were  calculated  from  dial  gauge 
readings  of  vertical  displacements.  Normal  loads  were  measured  by  both  a 
load  cell  and  an  air  pressure  gauge  on  the  ram.  Normal  load  was  compared  to 
the  sum  of  the  forces  through  the  upper  contacts  in  the  top  row  of  disks 
found  by  de  Josselin  de  Jong  and  Verruijt’s  photoelastic  force  measurement 
technique. 
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3.4  Image  Processing  System 

Both  still  and  moving  photoelastic  images  were  captured  and  processed 
by  a  digital  image  processing  and  analysis  system  consisting  of  the  following 
components:  (1)  a  Canon  LI  camcorder  with  a  410,000  pixel,  0.5-in.  CCD 
image  sensor  and  a  15  times  magnification  zoom  lens;  (2)  a  Macintosh  Quadra 
840 A V  desktop  computer  with  32MB  RAM,  1MB  VRAM  (video),  and  a 
Motorola  MC68040,  40  MHz  processor;  (3)  a  Data  Translation  DT2255  frame 
grabber  card;  and  (4)  NIH-Itnage^  image  processing  and  analysis  software. 

Both  qualitative  and  quantitative  analyses  were  performed  on 
photoelastic  images.  Full-field  images  of  the  entire  assembly  in  the  diffused- 
light  polariscope  were  qualitatively  analyzed  for  general  stress  transmission 
patterns  under  varying  degrees  of  simple  shear.  Individual  disk  contact  force 
magnitudes  and  directions  were  quantified  by  applying  the  modified  version 
of  de  Josselin  de  Jong  and  Verruijt's  technique  to  high-resolution  images 
from  the  lens  polariscope.  This  quantitative  analysis  was  used  to  produce 
force  vector  plots  for  comparison  with  results  from  simple  shear  experiments 
and  numerical  algorithms. 


1  NIH  Image  is  a  public  domain  program  and  was  written  by  Wayne  Rasband  at  the  U.S. 
National  Institute  of  Health.  It  is  available  from  the  Internet  by  anonymous  ftp  from 
zippy.nimh.nih.gov  or  on  floppy  disk  from  NTIS,  5285  Port  Royal  Rd.,  Springfield,  VA  22161, 
part  number  PB93-504868. 
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CHAPTER  4 

DEVELOPMENT  OF  CONTACT  FORCE  ANALYSIS  METHOD 

4.0  Introduction 

Paikowsky  and  DiRocco  (1993)  raised  some  concerns  regarding  the 
validity  of  the  photoelastic  method  of  contact  force  determination  proposed 
by  de  Josselin  de  Jong  and  Verruijt  (1969).  In  the  present  study  it  is  shown 
that  de  Josselin  de  Jong  and  Verruijt  (1969)  method  is  valid  provided  that 
several  adjustments  are  made  to  the  analysis  procedure  to  account  for  their 
assumptions.  In  this  chapter  the  assumptions  made  by  de  Josselin  de  Jong 
and  Verruijt  (1969)  are  discussed  and  a  modified  version  of  their  analysis 
technique  is  presented. 

4.1  De  Josselin  de  Jong  and  Verruijt  Method  of  Force  Analysis 

De  Josselin  de  Jong  and  Verruijt  (1969)  developed  a  relatively  simple 
technique  for  using  photoelastic  fringe  patterns  to  determine  the  magnitudes 
and  directions  of  contact  forces  between  cylinders.  They  made  use  of  the 
linear  relationship  between  the  contact  force  P  and  the  fringe  order  N  to 
determine  that 

p 

—  =  constant  (4.1) 

Nd 

where  d  is  the  diameter  of  the  circular  fringe  of  order  N.  This  constant  can  be 
calibrated  by  measuring  the  diameter  in  inches  of  a  circular  fringe  of  known 
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order  N  near  a  contact  point  where  a  known  force  P  acts.  This  value  is 
constant  for  any  model  made  from  the  same  photosensitive  material.  Once 
the  constant  is  known,  any  contact  force  can  be  determined  by  again  simply 
measuring  the  diameter  in  inches  of  a  circular  fringe  whose  order  is  known, 
and  solving  for  the  magnitude  of  P  from  Eq.  4.1. 

To  arrive  at  this  simple  constant  relationship,  however,  de  Josselin  de 
Jong  and  Verruijt  had  to  make  some  assumptions,  not  all  of  which  hold  true 
all  of  the  time.  They  first  assumed  that  the  fringes  that  develop  in  a  disk  are 
circles  which  pass  through  the  point  of  contact.  This  assumption  stems  from 
Frocht's  (1948)  work  on  photoelastic  analysis  of  a  concentrated  load  on  a  semi¬ 
infinite  plane,  that  is,  a  plane  which  extends  to  infinity  in  all  directions  on 
one  side  of  an  infinite  line  segment.  In  this  case,  the  fringes  that  develop  are 
indeed  exact  circles  that  pass  through  the  contact  point  (Fig.  4.1)  regardless  of 
the  angle  of  inclination  of  the  force. 


Figure  4.1.  Comparison  of  theoretical  and 
photoelastic  stress  patterns  (from  Frocht,  1948) 


A  contact  force  between  two  disks,  however,  is  not  exactly  represented 
by  a  concentrated  load  on  a  semi-infinite  plane,  which  can  therefore  be  only 
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an  approximation.  There  are  three  major  factors  that  separate  Frocht's  work 
with  a  concentrated  load  on  a  semi-infinite  plane  and  the  case  of  contact 
between  two  disks:  (1)  the  load  characteristics  (distributed  vs.  concentrated); 
(2)  the  plane  boundaries;  and  (3)  the  potential  superimposition  of  fringes 
from  nearby  forces  (Table  4.1). 

Table  4.1.  Differences  Between  the  Cases  of  (a)  a  Concentrated  Load  on  a  Semi-infinite  Plane 


and  (b)  Contact  Between  Two  Disks. 


Category 

(a)  Concentrated  Load  on  a 
Semi-Infinite  Plane 

(b)  Contact  Force  Between 
Two  Disks 

Load  Characteristics 

Concentrated 

Distributed 

Plane  Boundaries 

Infinite  Line  Segment 

Curved 

Potential  Superimposition  of 
Fringes  from  Nearby  Forces 

No 

Yes 

As  noted  by  Paikowsky  and  DiRocco  (1993),  these  differences  lead  to 
quite  different  fringe  geometries,  making  de  Josselin  de  Jong  and  Verruijt's 
assumptions  invalid.  However,  since  each  difference  affects  only  certain 
fringe  development  areas,  de  Josselin  de  Jong  and  Verruijt's  technique  can 
still  be  used  if  these  areas  are  avoided  when  choosing  a  fringe  for  analysis. 
Each  of  these  categories  are  discussed  below. 

4.1.1  Distributed  vs.  Concentrated  Loads 

Fringes  that  develop  as  a  result  of  concentrated  loads  have  significant 
geometric  differences,  in  the  region  near  the  point  of  contact,  from  those 
induced  by  distributed  loads  (Fig.  4.2). 
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Figure  4.2.  Close-up  photoelastic  images  of  disk  contacts,  (a) 


Concentrated  Load  (c)  Distributed  Load  (b,d)  Illustration  of  St.  Venant's 


Principle  (from  Frocht,  1948) 


Those  very  near  the  concentrated  load  are  so  close  together  that  they  begin  to 
blend  into  one  another  until  they  are  obscured  by  total  light  extinction.  Those 
very  near  the  distributed  load,  however,  are  quite  clear,  with  a  small. 
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ovalesque  fringe  surrounding  a  point  of  isotropy  near  the  contact.  The 
differences  in  fringe  geometry  that  result  from  the  two  different  loading  cases 
imply  that  the  choice  of  the  fringe  to  be  measured  in  de  Josselin  de  Jong  and 
Verruijt  technique  will  affect  the  calculated  magnitude  of  the  contact  force. 
To  resolve  this  issue,  St.  Venant's  principle  can  be  invoked. 

St.  Venant's  principle  states  that  stresses  produced  in  a  body  are 
identical  when  applied  by  any  statically  equivalent  load  system  except  over 
regions  near  the  points  of  load  application.  Therefore,  in  regions  sufficiently 
far  away  from  the  contact  point,  the  stresses  produced  by  the  distributed 
contact  load  are  identical  to  those  produced  by  a  concentrated  load. 

To  quantify  St.  Venant's  principle  for  the  specific  case  of  a  disk  in 
diametrical  compression,  Frocht  tested  the  photoelastic  disk  in  Figure  4.3  (the 
contacts  of  which  are  enlarged  in  Fig.  4.2)  with  both  large  (bottom)  and  small 
(top)  contact  areas  which  represent  approximately  distributed  and 
concentrated  loads,  respectively.  He  found  that  stresses  become  identical  at 
points  further  than  0.1  times  the  diameter  of  the  disk  from  the  contact  point. 


39 


Figure  4.3.  Disk  in  compression  by  pin  (top)  and  flat  surface 
(bottom)  to  approximate  concentrated  and  distributed  loads, 
(from  Frocht) 


De  Josselin  de  Jong  and  Verruijt's  technique  assumes  that  all  fringes 
are  circles  passing  through  the  point  of  contact.  From  Frocht's  tests  it  is  clear 
that  only  a  concentrated  load  produces  exactly  circular  fringes  in  the  region 
near  the  application  of  the  load.  But  by  St.  Venant's  principle  it  is  clear  that 
any  fringe  produced  by  a  disk  contact  (distributed)  force  that  is  at  least  O.ld 
from  the  contact  point  will  be  the  same  (circular)  fringe  that  would  have  been 
developed  by  the  application  of  a  concentrated  force.  Hence,  the  fringe  chosen 
for  analysis  must  be  at  least  O.ld  away  from  any  contact  to  be  free  from 
loading-induced  geometric  distortions. 
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4.1.2  Semi-infinite  Plane  vs.  Curving  Plane 

The  fringes  developed  in  a  semi-infinite  plane  from  a  concentrated 
force  will  remain  circular  at  any  distance  from  the  contact  point  (Fig.  4.1).  In  a 
disk,  however,  the  fringes  begin  to  be  affected  by  the  curving  boundary  of  the 
disk  (where  the  fringe  order  is  forced  to  remain  at  zero)  so  that  they  become 
oblong  rather  than  circular  (Fig.  4.4).  Again,  because  de  Josselin  de  Jong  and 
Verruijt's  technique  assumes  that  fringes  are  circular,  a  fringe  so  far  from  the 
contact  that  it  loses  its  circularity  cannot  be  used  in  analysis. 


Figure  4.4.  Disk  in  diametrical  compression  showing  diameter  d 
and  O.ld  zone. 


To  compare  fringes  developed  in  a  semi-infinite  plane  and  a  disk 
under  similar  loading  conditions,  a  disk  and  semi-infinite  plane  were 
compressed  together  (Fig.  4.5).  All  fringes  in  the  semi-infinite  plane  remain 
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circular  regardless  of  the  distance  from  the  contact,  while  the  disk's  fringes 


with  orders  smaller  than  about  3.5  become  distorted. 


Figure  4.5.  Disk  in  contact  with  half-plane  to  illustrate 
differences  in  fringe  development  between  a  disk  and  a  semi¬ 
infinite  plane.  Note  that  the  3.5  order  fringe  in  the  disk  is 
still  a  circle  and  is  outside  the  O.ld  zone. 


In  summary,  if  a  fringe  developed  by  two  disks  in  contact  is  too  near 
the  contact,  the  circular  geometry  that  would  have  resulted  from  a 
concentrated  load  will  be  distorted  by  the  distributed  load.  But  if  a  fringe  is 
too  far  from  the  contact,  its  circular  geometry  will  be  distorted  by  the  disk 
boundaries  so  that  it  becomes  an  ellipsoid.  The  fringe  chosen  for  analysis 
must  then  be  between  these  two  regions  if  it  is  to  remain  circular  as  presumed 
by  the  de  Josselin  de  Jong  and  Verruijt  technique. 
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4.1.3  Superimposition  of  Fringes 


Paikowsky  and  DiRocco  (1993)  addressed,  in  considerable  detail,  the 
concern  that  since  the  fringe  pattern  is  a  result  of  the  superimposition  of 
stresses  induced  by  all  forces  acting  on  the  disk,  the  circular  shape  of  a  fringe 
induced  by  one  force  will  be  altered.  This  is  true  in  regions  further  away  from 
contact  points,  where  multiple  forces  influence  fringe  shape,  but  the 
influence  of  a  force  in  the  region  near  its  contact  point  dramatically 
overwhelms  the  influence  of  other  relatively  near  forces. 

Paikowsky  and  DiRocco  performed  tests  to  illustrate  fringe 
superimposition.  They  loaded  a  1.2"  disk  with  two  different  loading  schemes 
(Fig.  4.6a).  The  disk  used  was  1.2"  in  diameter,  which  means  that  the  fringe 
chosen  for  analysis  by  de  Josselin  de  Jong  and  Verruijt's  technique,  as 
superimposition  becomes  a  factor  further  away  from  the  contact,  but  does  not 
affect  fringes  in  the  area  of  interest. 


Figure  4.6.  (a)  Loading  configurations  for  1.2"  disk  with  both  three  loads  and 
diametrical  loads,  (b)  Enlargement  of  fringes  developed  near  upper  right 
contact  for  both  cases.  Solid  line  represents  three  contact  loading,  broken  line 
represents  diametrical  loading.  Fringes  are  almost  identical  in  the  region 
O.ld  (.12")  from  contact,  (from  Paikowsky  and  DiRocco,  1993) 


44 


The  loading  situations  that  are  investigated  for  the  current  work 
contain  small  enough  loads  that  neighboring  forces  are  not  expected  to  have 
significant  influence  upon  the  fringes  in  the  regions  used  for  analysis. 

It  is  concluded  that  de  Josselin  de  Jong  and  Verruijt  technique  for 
determining  contact  forces  simply  by  geometric  measurement  of  fringes  is 
indeed  valid,  provided  that  only  fringes  within  a  certain  region  are 
considered  in  analysis.  Fringes  measured  for  contact  force  determination 
must  have  a  diameter  of  approximately  0.1  times  the  diameter  of  the  disk  and 
must  pass  through  the  contact  point  under  analysis. 
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4.2  Application  of  Force  Analysis  Method 

From  a  photoelastic  image  of  two  right  circular  cylinders  (disks)  in 
contact,  it  is  possible  to  determine  both  the  magnitude  and  direction  of  the 
contact  force  between  them,  provided  that  model  material  has  been  calibrated. 
By  slightly  modifying  the  method  proposed  by  de  Josselin  de  Jong  and 
Verruijt  (1971)  for  force  determination  through  a  single  contact,  and  using  an 
image  processing  system  to  automate  several  time-consuming  steps,  the  load 
transmission  pattern  through  a  photoelastic  model  of  a  granular  media  can  be 
determined. 

Using  a  camcorder  to  capture  digital  images  and  a  framegrabber  card  to 
store  them  on  a  computer  eliminates  the  former  need  to  shoot  and  develop 
careful  photographs  for  later  hand  analysis.  The  image  processing  software 
allows  geometric  measurements  such  as  distances  and  angles  to  be  made 
quickly  and  accurately. 

After  the  disks  are  packed  into  the  loading  device  and  subjected  to  load 
and  shear,  the  assembly  is  placed  in  line  with  the  polariscope.  The  camera  is 
positioned  to  capture  high-resolution  close-ups  of  each  disk.  All  images, 
including  a  full-field  shot  and  close-ups,  are  captured  at  one  sitting  and  stored 
for  later  analysis.  After  all  images  have  been  captured,  the  force  magnitude 
and  direction  at  each  contact  within  the  assembly  is  obtained  by  the  following 
procedure. 
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4.2.1  Contact  Force  Magnitude 

The  full-field  image  (Fig.  4.7)  is  captured  while  the  assembly  is  in  the 
diffused-light  polariscope  for  its  greater  viewing  area.  This  image  gives  the 
general,  qualitative  stress  transmission  pattern  and  is  used  for  later 
comparison  with  a  force  vector  plot  made  from  the  quantification  of  forces  at 
individual  contacts. 


Figure  4.7.  Full-field  image  of  a  dense  packing  of  0.75-in.  diameter  disks  under  145  lb.  normal 
load  and  a  shear  strain  of  0.0186.  Captured  from  a  diffused-light  polariscope. 


Individual  disk  images  are  then  taken  from  the  lens  polariscope  for  its 
sharper,  more  detailed  images  (Fig.  4.8).  NIH  Image  1.57  (Image),  the  image 
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processing  software  used  for  the  current  work,  is  then  used  to  manipulate  the 
image  for  all  necessary  measurements.  The  diameter  of  the  disk,  known  in 
this  case  to  be  0.75  in.,  is  measured  in  pixels,  and  the  scale  of  the  image  is 
calculated  and  recorded. 


Figure  4.8.  Close-up  image  of  a  0.75-m.  diameter  disk  from  dense  arrangement  in  Fig.  4.7 
(bottom  row,  fifth  from  left).  Note  the  improved  contrast  (darker  fringes)  in  this  image 
captured  from  a  lens  polariscope. 


The  contact  to  be  measured  is  magnified  for  better  manipulation.  A 
circular  region  of  interest  (ROI)  equal  to  0.1  times  the  diameter  of  the  disk  is 
created  and  placed  so  that  it  approximates  a  contact  force-induced  fringe  (Fig. 
4.9).  (It  was  shown  above  that  a  fringe  produced  by  a  distributed  force  on  a 
disk  boundary  will  have  the  same  circular  geometry  as  a  fringe  produced  by  a 


48 


concentrated  force  on  a  semi-infinite  plane  provided  that  the  fringe  diameter 
is  approximately  0.1  times  the  diameter  of  the  disk,  and  passes  through  the 
contact  in  question.)  If  a  fringe  does  not  fall  exactly  on  the  ROI,  either  (a)  the 
fringe  order  at  the  ROI  can  be  interpolated  from  the  two  bordering  fringes,  or 
(b)  the  ROI  can  be  slightly  reduced  or  enlarged  to  meet  the  nearest  fringe. 


Figure  4.9.  Magnification  of  upper  left  contact  point  of  disk  in  Figure  2. 
ROI  with  diameter  of  0.075,  or  O.ld,  is  shown  to  coincide  with  fringe  of 
order  1.5.  Force  can  be  calculated  from  the  equation:  P/Nd  =  constant. 


It  was  shown  that  P/Nd  is  a  constant  for  any  given  material  (Eq.  4.1), 
where  P  is  the  force  through  the  contact,  N  is  the  fringe  order,  and  d  is  the 
diameter  of  the  circle  that  best  approximates  the  fringe.  After  this  constant 
has  been  calibrated  for  the  model  material,  the  force  through  the  fringe  is 
easily  calculated  from  the  constant. 
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4.2.2  Contact  Force  Direction 


It  is  also  known  that  a  contact  force  passes  through  the  centroid  of  any 
contact-induced  fringe  whose  diameter  is  greater  than  one-tenth  that  of  the 
disk.  Within  the  O.ld  region,  fringes  are  often  distorted  by  non-uniform 
boundary  conditions  that  are  intensified  by  distributed  loading,  and  should  be 
avoided.  The  centroid  of  a  circular  fringe  that  is  approximated  by  an  ROI  is 
easily  calculated  by  Image.  The  line  segment  that  connects  the  contact  point 
and  the  centroid  of  the  measured  fringe  coincides  with  the  line  of  action  of 
the  force  (Fig.  4.10). 


Figure  4.10.  Centroid  of  fringe  marked  with  a  plus  and  contact 
force  vector  plotted  with  origin  at  the  contact  and  passing 
through  fringe  centroid.  Angle  with  contact  plane  can  easily  be 
calculated  with  Image. 
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4.3  Semi-automation 

The  above  procedure  has  been  semi-automated  by  programs  written  to 
take  advantage  of  several  of  Image’s  built-in  measurement  features. 
Routines  were  written  that  allow  incremental  manipulation  of  the  size  and 
position  of  the  ROI,  that  automatically  calculate  the  force  with  manual  input 
of  only  the  fringe  order,  and  that  find  and  mark  the  center  of  the  ROI  for 
determination  of  contact  force  direction.  If  large  assemblies  with  many 
contacts  are  to  be  analyzed,  it  may  be  worthwhile  to  attempt  the  automation 
of  fringe  order  determination  based  on  coordinates  of  the  disk  boundary  and 
contacts,  and  fringe/backgroimd  brightness  contrast.  The  assemblies  analyzed 
for  this  work  were  small  enough  that  manual  analysis  was  less  time- 
consuming  than  attempting  to  program  detailed  automation.  The  semi¬ 
automation,  however,  allowed  more  rapid  measurement  and  recording  than 
would  have  been  achieved  by  hand-measuring  photographic  prints. 

The  force  magnitude  and  direction  for  each  fringe  is  then  recorded  in  a 
file,  and  used  to  plot  the  force  vector  diagram  (Chapter  5)  for  the  assembly. 
Creation  of  a  force  vector  diagram  is  also  a  procedure  that  could  potentially  be 
automated  within  Image  if  desired,  but  was  not  done  for  this  work  for  the 
reasons  mentioned  above. 
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CHAPTER  5 

RESULTS  AND  DISCUSSION 

5.0  Introduction 

Several  qualitative  observations  have  been  made  regarding  the 
behavior  of  the  disk  assemblies  in  simple  shear.  These  include  the  frictional 
behavior  between  disks,  disks  rotations,  volumetric  strains,  boundary  stress 
uniformity,  and  load  distribution.  Force  vector  plots  have  been  obtained  for 
different  stages  of  simple  shear  deformation  of  the  disk  assemblies. 

5.1  Behavioral  Observations 

5.1.1  Inter  disk  Frictional  Behavior 

Motion  pictures  reveal  that  a  slipping  and  locking  mechanism  exists  at 
disk  contacts.  As  the  disks  are  forced  to  move  relative  to  one  another,  their 
surfaces  either  slide  past  or  lock  together,  the  latter  occurring  when  the 
interface  friction  between  the  disk  surfaces  is  sufficient  to  resist  the  force 
driving  the  relative  movement. 

Throughout  the  range  of  the  shear  deformation,  the  disks 
continuously  engage  and  disengage  one  another.  Interdisk  forces  lock  the 
disks  together,  then  overcome  the  friction  between  them  by  either  growing 
too  large  for  the  friction  to  resist,  or  by  acting  at  an  angle  larger  than  the 
maximum  angle  of  friction  resistance.  In  a  large  granular  assembly  such  as 
sand,  this  mechanism  is  quantified  by  assigning  the  material  an  angle  of 


internal  friction,  an  average  value  that  is  affected  by  grain  shape,  size, 
distribution  and  constituent  minerals. 

5.1.2  Volumetric  Strain 

Granular  materials  have  been  observed  to  undergo  volume  change 
when  sheared  in  triaxial,  direct  shear,  and  simple  shear  tests.  Densely  packed 
materials  tend  to  dilate  (increase  in  volume),  while  loosely  packed  materials 
tend  to  compress  (decrease  in  volume).  Dense  packings  exhibit  significant 
particle  interlocking  which  must  be  overcome  when  sheared,  forcing  the 
particles  to  roll  up  and  over  one  another,  increasing  the  volume  of  the 
specimen.  Loose  packings  show  less  particle  interlocking,  so  the  grains  can 
slide  into  void  spaces,  resulting  in  compression. 

In  the  present  study,  extremes  of  both  conditions  are  examined.  The 
loose  assembly  of  disks,  described  as  a  simple  cubic  packing,  represents  the 
loosest  way  in  which  the  disks  can  be  arranged,  one  stacked  directly  upon 
another.  The  dense  assembly,  described  as  a  simple  stagger,  is  the  densest 
possible  arrangement  of  disks  (Fig.  5.1). _ 


Figure  5.1.  Photoelastic  disk  arrangements.  Simple  cubic  and  simple  stagger  packings 
representing  loose  (eo  =  0.27)  and  dense  (eo  =  0.15)  granular  assemblies,  respectively. 
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Volumetric  strain  is  often  plotted  against  axial  strain  in  triaxial  and 
direct  shear  tests  or  against  shear  strain  in  simple  shear  tests.  Volumetric 
strain,  Sy,  is  defined  as 


£v  = 


(5.1) 


where  AV  is  the  change  in  volume  of  the  specimen  and  Vo  is  its  initial 
volume.  Since  only  vertical  movement  of  assembly  boundaries  is  allowed, 
AV  =  Ay  where  Ay  is  the  vertical  movement  of  the  top  platen,  and  Vo  =  ho 
where  ho  is  the  initial  height  of  the  assembly.  Volumetric  strain  versus  shear 
strain  plots  for  both  the  dense  and  loose  arrangements  of  disks  follows 
similar  patterns  to  those  of  sands  (Fig.  5.2). 
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Volumetric  Strain  vs.  Shear  Strain 


Figure  5.2.  Volumetric  strain  vs.  shear  strain.  Comparison  of  volumetric 
changes  with  shear  strain  for  dense  and  loose  disk  arrangements. 

Normally,  a  dense  sand  will  exhibit  some  compression  during  the 
initial  stages  of  shearing,  and  will  strongly  dilate  thereafter.  Since  the  dense 
disk  assembly  is  in  its  densest  possible  state,  it  shows  no  initial  compression 
and  exhibits  immediate  dilatation.  Likewise,  a  loose  sand,  after  some 
shearing,  will  become  denser  by  compression  and,  because  of  its  increased 
density,  will  begin  to  show  a  lower  rate  of  compression.  The  loose  disk 
assembly  is  in  its  loosest  state,  however,  and  therefore  requires  significant 
shearing  to  reach  a  dense  enough  state  to  begin  to  exhibit  the  reduced 
compression  rate. 
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5.1.3  Disks  Rotations 

Disks  rotations  in  a  particulate  medium,  depending  on  the  degree  of 
particle  angularity,  can  significantly  affect  the  path  by  which  a  load  is 
transferred  through  the  medium.  Understanding  why  and  how  individual 
disk  rotates  during  simple  shear  is  an  integral  part  of  an  attempt  to  model  its 
behavior. 

To  measure  the  rotation  of  individual  disks  during  simple  shear,  a 
reference  line  (appearing  as  a  dark  line  in  the  photoelastic  images)  was  etched 
across  the  diameter  of  a  representative  disk  in  the  assembly,  and  its  position 
was  measured  before  and  after  a  desired  simple  shear  strain  was  applied  (Figs. 
5.3, 5.4, 5.5). 
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Figure  5.3.  Rotation  of  grains,  (a)  Dark  line  on  disk  in  pre-shear  reference  position,  (b)  Line 
in  post-shear  position  indicates  11°  clockwise  rotation  during  shear.  Photoelastic  image 
shows  no  stress  band  development  on  the  vertical  planes  indicating  that  no  locking  with 
neighboring  disks  is  occurring  on  these  planes. 
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In  the  loose  packing,  which  tends  to  deform  as  rigid  columns  of 
particles  sliding  relative  to  one  another,  individual  disk  rotation  was 
observed  in  some  cases.  In  general,  if  the  vertical  contacts  of  a  disk  in  this 
arrangement  was  sliding  relative  to  neighboring  disks,  the  individual  particle 
rotated  relative  to  the  loading  frame  (Fig.  5.3).  If  the  vertical  contacts  are 
exhibiting  friction  locking,  the  disk  will  exhibit  much  less  rotation  relative  to 
the  frame  (Fig.  5.4). 


Figure  5.4.  Rotation  of  grams,  (a)  Dark  line  on  disk  in  pre-shear  reference  position,  (b)  Line 
in  post-shear  position  indicates  negligible  rotation.  The  stress  band  development  on  the 
vertical  contacts  indicates  locking  with  these  neighboring  disks.  This  mechanism  is 
hypothesized  to  be  responsible  for  the  difference  in  rotation  behavior  from  the  disk  in  Fig. 
5.3. 


A  disk  in  the  densest  possible  arrangement  has  been  observed  to 
remain  without  rotating  relative  to  the  frame  regardless  of  whether  it  is 
sliding  or  locking  with  neighboring  disks  (Fig.  5.5). 


positions,  undergoes  no  rotation  during  shear. 


5.1.4  Load  Distribution 

The  observed  pattern  of  stress  transmission  in  the  loose  disk 
arrangement  under  simple  shear  is  somewhat  random.  Before  shearing,  the 
disks  are  arranged  in  eight  columns  of  six  disks  each,  and  the  vertical  load 
applied  by  the  top  platen  is  transmitted  vertically  through  each  column  (Fig. 
5.6). 

As  shear  is  applied,  the  disks  begin  to  slide  relative  to  one  another  at 
those  contacts  where  the  interdisk  friction  is  overcome,  and  lock  together  at 
those  where  the  friction  provided  between  disks  is  enough  to  resist  the 
driving  contact  force.  The  selection  of  these  contacts  appears,  from  resulting 
photoelastic  images,  to  be  quite  arbitrary.  Three  tests  ran  under  the  identical 
vertical  loads  and  shear  strains  exhibited  significantly  different  stress 
transmission  patterns  (Figs.  5.7a-c). 


Figure  5.6.  Loose  packing  before  shear  strain  application.  Notice  that  load  is  transferred 
vertically  through  each  column  of  disks. 


Figure  5.7a.  Photoelastic  image  of  loose  assembly  under  simple  shear.  Note  that  shear 
changes  the  stress  transmission  pattern  nonuniformly  such  that  some  vertically  oriented 
disk  contacts  show  load  transfer  while  others  do  not. 


Figure  5.7b.  A  second  photoelastic  image  of  loose  assembly  imder  simple  shear.  Test 
procedure  was  exactly  the  same  as  in  Fig.  5.7a,  applying  identical  vertical  loads  and  shear 
strains.  Note  that  the  stress  is  transferred  across  different  vertical  contacts  in  the  two 
figures. 


Figure  5.7c.  A  third  photoelastic  image  of  loose  assembly  imder  the  same  vertical  load  and 
shear  strain  as  in  Figs.  5.7a  and  5.7b.  Stress  transmission  pattern  is  again  transferred  across 
different  vertical  contacts. 

This  is  partially  attributable  to  the  nature  of  the  loose  packing.  In 
simple  shear,  the  vertical  load  that  initially  caused  stress  to  flow  vertically  is 
redirected  to  flow  across  the  box's  shorter  diagonal.  In  the  idealized  loose 
packing,  the  disks  are  arranged  such  that  the  contacts  are  90°  from  one 
another,  making  flow  impossible  between  two  diagonal  (i.e.,  non-contacting) 
disks.  The  expected  diagonal  flow  pattern  is,  however,  exhibited  by  the  dense 
packing.  Load  distribution  in  the  dense  assembly  is  discussed  in  detail  in 


Section  5.2.1. 
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Drescher  and  de  Josselin  de  Jong  (1972)  explained  this  phenomenon 
with  a  "double-sliding  free-rotating  model",  in  which  the  two  mechanisms  of 
rigid-block  sliding  and  rotation  are  separate  but  simultaneous.  This  model 
was  derived  from  Horne's  (1965)  observation  that  the  deformation  of  a 
granular  material  results  from  the  relative  movement  of  instantaneously 
rigid  groups  of  particles,  that  constantly  redistribute  themselves  in  division 
and  coalescence. 

There  are  infinite  combinations  in  which  the  assembly  of  disks  can 
form  rigid,  sliding  blocks  (Fig  5.8a,b),  then  rotate,  if  necessary,  to  retain  the 
geometric  boundary  constraints  of  the  assembly.  From  all  combinations,  the 
one  that  provides  the  least  resistance  is  chosen. 


Sliding 

(a) _ 

Figure  5.8.  Bookstack  analogy  of  simple  shear 
deformation  (from  de  Josselin  de  Jong  (1972). 


Rotation 


(b) 


Microscopic  variations  on  the  surface  of  the  disks,  left  behind  from 
machining,  produce  a  non-uniform  roughness  over  the  surface  of  one  disk. 
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and  a  non-constant  roughness  from  disk  to  disk.  These  variations  provide 
the  frictional  variation  throughout  the  assembly  required  to  produce 
inconsistent  deformation  combinations  from  test  to  test  in  an  apparently 
homogeneous  assembly.  Such  microscopic  inconsistencies  are  not  accounted 
for  in  continuum  mechanics  but  which  can  significantly  affect  load 
distribution  calculations. 

5.1.5  Uniformity  of  Boundary  Stress 

Theoretical  approaches  based  on  elasticity  (Roscoe,  1953)  and 
elastoplasticity  (Budhu  and  Britto,  1987),  and  experimental  approaches 
(Stroud,  1971;  Budhu,  1979)  have  shown  that  nonuniformities  develop  along 
boundaries  of  specimens  undergoing  simple  shear  deformation. 
Nonuniform  distribution  of  both  shear  and  normal  stresses  along  the  top  and 
bottom  boundaries  of  sand  specimens  was  measured  in  the  Cambridge 
University  simple  shear  device  which  was  equipped  with  contact  stress 
transducers,  (Budhu,  1983). 

Photoelastic  analysis  can  detect  nonuniform  stress  development. 
Images  of  the  dense  assembly  of  photoelastic  disks  showed  that  the  normal 
stress  nonuniformity  is  similar  to  that  observed  in  tests  on  sand  (Fig.  5.9). 
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Figure  5.9.  Boundary  nonuniformities.  Left:  Photoelastic  image  of  a  dense  granular  assembly 
rmder  simple  shear.  Note  that  disks  at  top  right  and  bottom  left  show  no  stress,  while  the 
opposite  comers  exhibit  high  stresses.  Right:  Diagram  showing  normal  stress  nonxmiformities 
of  a  simple  shear  test  on  sand  (from  Budhu,  1980).  Experiments  on  disks  and  on  sand 
qualitatively  agree. 

As  the  specimen  is  sheared,  the  stress  flow,  which  was  initially  vertically 
oriented,  gradually  shifts  to  one  that  flows  along  the  shortest  diagonal  of  the 
shear  box.  This  is  reflected  in  stress  distributions  on  the  top  and  bottom 
boundaries  wherein  higher  stresses  developed  on  the  top  left  and  bottom 
right  comers. 
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5.2  Comparison  With  the  Lattice  Model 

On  application  of  the  vertical  load  with  no  applied  shear  strain,  the 
internal  distribution  of  stresses  in  the  dense  disk  assembly,  as  revealed  by  the 
photoelastic  image  (Fig.  5.10),  was  not  uniform.  Under  a  relatively  small 
shear  strain  (but  large  enough  to  exhibit  visible  directional  stress  flow)  the 
disks  formed  rigid  diagonal  columns  that  transmit  the  applied  normal  force 


Figure  5.10.  Photoelastic  image  of  dense  packing  with  no  shear  strain. 


Figure  5.11.  Photoelastic  image  of  dense  packing  with  shear  strain  (0.1%). 

Force  vector  plots  (Fig.  5.12a)  from  the  photoelastic  image  (Fig.  5.10)  showed 
that  the  vertical  applied  force  is  mainly  transmitted  by  two  pairs  of  inclined 
columns  of  disks  originating  from  the  two  extreme  pair  of  disks  at  the  top 
right  and  left  corners  and  meeting  at  the  middle  pair  of  disks  at  the  bottom 
boundary.  The  lattice  model  predicted  a  similar  force  distribution  (Fig.  5.12b) 
except  that  rather  than  two  pairs  of  columns  there  are  two  distinct  columns 
originating  from  the  disks  at  the  extreme  corners  of  the  top  boundary  meeting 
at  the  middle  pair  of  disks  at  the  bottom  boundary.  These  results  are 
intuitively  correct  since  the  disks  at  the  extreme  comers  of  the  top  boundary 
have  only  two  neighboring  disks  each  to  transmit  the  vertical  load  compared 
to  four  disks  for  the  other  disks  on  the  top  boundary. 
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Figure  5.12a.  Contact  force  vector  plot  of  dense  assembly  at  zero  shear  strain  (Fig 
5.10)  created  by  using  the  modified  version  of  de  Josselin  de  Jong  and  Verruijt's 
method  of  finding  contact  force  vectors. 


Figure  5.12b.  Contact  force  vector  plot  at  zero  applied  shear  strain  (Fig  5.10) 
predicted  by  the  lattice  type  model. 
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The  nonuniform  distribution  of  load  is  bound  to  occur  from  the  structural 
arrangement  of  the  media  and  the  conditions  at  the  boundaries.  One  cannot 
expect  to  obtain  exact  replica  of  the  force  distribution  from  the  simple  shear 
tests  and  the  lattice  model.  Although  the  disks  used  in  the  tests  were 
machined  with  close  tolerances,  even  small  surface  variations  can 
significantly  alter  the  force  distribution  from  an  ideal  array  of  disks  as 
assumed  in  the  LTM.  As  the  shear  strain  is  applied  to  the  assembly,  load 
redistribution  occurs  and  again  LTM  predicts  load  paths  that  are  similar  to  the 
experiments  (Figs.  5.13a,b). 

For  the  loose  assembly,  the  disks  are  aligned  one  on  top  of  the  other  so 
that  the  load  paths  are  vertical  (Fig.  5.6)  under  the  vertical  load  only.  The 
LTM  results  show  a  similar  load  path  (Fig.  4.9b,  PART  I).  These  initially 
vertical  paths  (columns)  become  inclined  as  shear  strains  are  applied 
(compare  Fig.  5.7b  with  Fig.  4.9b,  PART  I).  Because  of  the  slight  discrepancies 
in  the  size  and  contact  roughness  of  the  disks  in  the  experiments,  it  is  not 
possible  to  conduct  a  realistic  comparison  between  the  experimental  results 
and  the  results  of  the  LTM.  However,  the  predictions  from  the  LTM  and  the 
experiments  are  qualitatively  similar. 
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Figure  5.13a.  Contact  force  vector  plot  dense  assembly  at  a  small  shear  strain 
(0.1%)  (Fig  5.11)  created  by  using  the  modified  version  of  de  Josselin  de  Jong  and 
Verruijt’s  method  of  finding  contact  force  vectors. 


Figure  5.13b.  Contact  force  vector  plot:  dense  assembly  at  a  small  shear  strain  (Fig 
5.11)  predicted  by  the  lattice  type  model. 


CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions 

An  experimental  system  consisting  of  a  shear  box,  photoelastic 
equipment  and  digital  image  processing  technology  was  developed  to 
investigate  load  distribution  and  preferred  deformation  mechanisms 
involved  in  the  deformation  of  a  particulate  media  under  simple  shear.  The 
system  was  very  effective.  Micro-features  that  affect  the  behavior  of 
particulate  media  can  be  individually  addressed  and  analyzed  in  great  detail. 

A  method  of  contact  force  analysis  based  on  photoelastic  data  has  been 
developed  and  implemented  to  determine  contact  force  vector  distribution 
networks  for  an  assembly  of  photoelastic  disks.  Comparison  of  the  load 
distribution  networks  with  similar  networks  predicted  by  the  lattice  type 
model  reveal  concordance.  The  lattice  type  model  appears  to  be  a  viable 
model  to  provide  insights  into  the  micromechanical  behavior  of  granular 
assemblies. 

Photoelastic  motion  pictures  reveal  that  a  slipping  and  locking 
mechanism  exists  at  disk  contacts.  Throughout  the  range  of  the  shear 
deformation,  the  disks  continuously  engage  and  disengage  one  another. 
Clusters  or  rigid  groups  of  particles  are  formed  during  deformation.  This  is 
consistent  with  Horne's  (1965)  observation  that  the  deformation  of  a  granular 
material  results  from  the  relative  movement  of  instantaneously  rigid  groups 
of  particles,  constantly  redistributing  themselves  in  division  and  coalescence. 
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It  is  also  consistent  with  Drescher  and  de  Josselin  de  Jong's  (1972)  double- 
sliding,  free-rotating  model. 

Rotation  of  individual  disk  during  simple  shear  has  been  observed  to 
be  dependent  upon  whether  the  disk  in  question  is  sliding  relative  to 
neighboring  particles.  In  general,  if  the  vertical  contacts  of  a  disk  in  the  loose 
arrangement  was  slipping  relative  to  neighboring  disks,  the  individual  disk 
rotated  relative  to  the  loading  frame.  If  the  vertical  contacts  exhibited  friction 
locking,  the  disk  showed  almost  no  rotation  relative  to  the  frame.  A  disk  in 
the  densest  possible  arrangement  has  been  observed  to  remain  without 
rotating  relative  to  the  frame,  probably  because  it  has  more  contact  points 
with  neighboring  disks,  increasing  the  likelihood  that  it  is  locked  against  one 
or  more  of  those  disks. 

The  pattern  of  stress  transmission  in  the  loose  disk  arrangement  under 
simple  shear  has  been  observed  to  be  somewhat  random.  This  is  partially 
attributable  to  the  nature  of  the  loose  packing.  In  simple  shear,  the  vertical 
load  that  initially  caused  stress  to  flow  vertically  is  redirected  to  flow  across 
the  box's  shorter  diagonal.  In  the  idealized  loose  packing,  the  disks  are 
arranged  such  that  the  contacts  are  90°  from  one  another,  making  load 
transfer  impossible  between  two  diagonal  (i.e.,  non-contacting)  disks.  Load 
redistribution  during  simple  shear  is  sensitive  to  frictional  micro- variations 
in  the  surface  profile  of  the  particles. 
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6.2  Recommendations 

1.  The  photoelastic  experiments  coupled  with  the  image  processing  software 
is  a  powerful  experimental  method  to  determine  load  paths  and  deformation 
mechanisms.  This  arrangement  can  be  used  to  make  more  direct 
comparisons  with  results  from  the  lattice  type  model  and  discrete  element 
methods  provided  that  the  disks  are  precisely  the  same  size  and  roughness. 
This  is  possible  by  molding  the  disks  rather  than  cutting  them  out  from 
photoelastic  sheets. 

2.  A  better  simulation  can  be  accomplished  using  photoelastic  spheres  in  a 
modified  version  of  the  simple  shear  device  described  here.  The 
modification  involves  the  redesign  of  the  shear  box  to  accommodate  the 
spheres. 
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